OpenCPN Partial API docs
Loading...
Searching...
No Matches
chartimg.cpp
1/**************************************************************************
2 * Copyright (C) 2015 by David S. Register *
3 * *
4 * This program is free software; you can redistribute it and/or modify *
5 * it under the terms of the GNU 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 General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program; if not, see <https://www.gnu.org/licenses/>. *
16 **************************************************************************/
17
24#include <assert.h>
25
26// For compilers that support precompilation, includes "wx.h".
27#include <wx/wxprec.h>
28
29#ifndef WX_PRECOMP
30#include <wx/wx.h>
31#endif
32
33// Why are these not in wx/prec.h?
34#include <wx/dir.h>
35#include <wx/stream.h>
36#include <wx/wfstream.h>
37#include <wx/tokenzr.h>
38#include <wx/filename.h>
39#include <wx/image.h>
40#include <wx/fileconf.h>
41
43
44#include "config.h"
45#include "chartimg.h"
46#include "ocpn_pixel.h"
47
48#ifndef _MSC_VER
49#include <signal.h>
50#include <setjmp.h>
51#include "navutil.h"
52
53#define OCPN_USE_CONFIG 1
54#endif
55
56#ifndef _MSC_VER
57struct sigaction sa_all_chart;
58struct sigaction sa_all_previous;
59
60sigjmp_buf env_chart; // the context saved by sigsetjmp();
61
62void catch_signals_chart(int signo) {
63 switch (signo) {
64 case SIGSEGV:
65 siglongjmp(env_chart, 1); // jump back to the setjmp() point
66 break;
67
68 default:
69 break;
70 }
71}
72
73#endif
74
75// Missing from MSW include files
76#ifdef _MSC_VER
77typedef __int32 int32_t;
78typedef unsigned __int32 uint32_t;
79typedef __int64 int64_t;
80typedef unsigned __int64 uint64_t;
81#endif
82
83// ----------------------------------------------------------------------------
84// Random Prototypes
85// ----------------------------------------------------------------------------
86
87bool G_FloatPtInPolygon(MyFlPoint *rgpts, int wnumpts, float x, float y);
88
89// ----------------------------------------------------------------------------
90// private classes
91// ----------------------------------------------------------------------------
92
93// ============================================================================
94// ThumbData implementation
95// ============================================================================
96
97ThumbData::ThumbData() { pDIBThumb = NULL; }
98
99ThumbData::~ThumbData() { delete pDIBThumb; }
100
101// ============================================================================
102// Palette implementation
103// ============================================================================
104opncpnPalette::opncpnPalette() {
105 // Index into palette is 1-based, so predefine the first entry as null
106 nFwd = 1;
107 nRev = 1;
108 FwdPalette = (int *)malloc(sizeof(int));
109 RevPalette = (int *)malloc(sizeof(int));
110 FwdPalette[0] = 0;
111 RevPalette[0] = 0;
112}
113
114opncpnPalette::~opncpnPalette() {
115 if (NULL != FwdPalette) free(FwdPalette);
116 if (NULL != RevPalette) free(RevPalette);
117}
118
119// ============================================================================
120// ChartBase implementation
121// ============================================================================
122ChartBase::ChartBase() {
123 m_depth_unit_id = DEPTH_UNIT_UNKNOWN;
124
125 pThumbData = new ThumbData;
126
127 m_global_color_scheme = GLOBAL_COLOR_SCHEME_RGB;
128
129 bReadyToRender = false;
130
131 Chart_Error_Factor = 0;
132
133 m_Chart_Scale = 10000; // a benign value
134 m_Chart_Skew = 0.0;
135
136 m_nCOVREntries = 0;
137 m_pCOVRTable = NULL;
138 m_pCOVRTablePoints = NULL;
139
140 m_nNoCOVREntries = 0;
141 m_pNoCOVRTable = NULL;
142 m_pNoCOVRTablePoints = NULL;
143
144 m_EdDate = wxInvalidDateTime;
145
146 m_lon_datum_adjust = 0.;
147 m_lat_datum_adjust = 0.;
148
149 m_projection = PROJECTION_MERCATOR; // default
150}
151
152ChartBase::~ChartBase() {
153 delete pThumbData;
154
155 // Free the COVR tables
156
157 for (unsigned int j = 0; j < (unsigned int)m_nCOVREntries; j++)
158 free(m_pCOVRTable[j]);
159
160 free(m_pCOVRTable);
161 free(m_pCOVRTablePoints);
162
163 // Free the No COVR tables
164
165 for (unsigned int j = 0; j < (unsigned int)m_nNoCOVREntries; j++)
166 free(m_pNoCOVRTable[j]);
167
168 free(m_pNoCOVRTable);
169 free(m_pNoCOVRTablePoints);
170}
171
172wxString ChartBase::GetHashKey() const {
173 wxString key = GetFullPath();
174 wxChar separator = wxFileName::GetPathSeparator();
175 for (unsigned int pos = 0; pos < key.size(); pos = key.find(separator, pos))
176 key.replace(pos, 1, "!");
177 return key;
178}
179
180/*
181int ChartBase::Continue_BackgroundHiDefRender()
182{
183 return BR_DONE_NOP; // signal "done, no refresh"
184}
185*/
186
187// ============================================================================
188// ChartDummy implementation
189// ============================================================================
190
191ChartDummy::ChartDummy() {
192 m_pBM = NULL;
193 m_ChartType = CHART_TYPE_DUMMY;
194 m_ChartFamily = CHART_FAMILY_UNKNOWN;
195 m_Chart_Scale = 22000000;
196
197 m_FullPath = "No Chart Available";
198 m_Description = m_FullPath;
199}
200
201ChartDummy::~ChartDummy() { delete m_pBM; }
202
203InitReturn ChartDummy::Init(const wxString &name, ChartInitFlag init_flags) {
204 return INIT_OK;
205}
206
207void ChartDummy::SetColorScheme(ColorScheme cs, bool bApplyImmediate) {}
208
209ThumbData *ChartDummy::GetThumbData(int tnx, int tny, float lat, float lon) {
210 return (ThumbData *)NULL;
211}
212
213bool ChartDummy::UpdateThumbData(double lat, double lon) { return FALSE; }
214
215bool ChartDummy::GetChartExtent(Extent *pext) {
216 pext->NLAT = 80;
217 pext->SLAT = -80;
218 pext->ELON = 179;
219 pext->WLON = -179;
220
221 return true;
222}
223
224bool ChartDummy::RenderRegionViewOnGL(const wxGLContext &glc,
225 const ViewPort &VPoint,
226 const OCPNRegion &RectRegion,
227 const LLRegion &Region) {
228 return true;
229}
230
231bool ChartDummy::RenderRegionViewOnDC(wxMemoryDC &dc, const ViewPort &VPoint,
232 const OCPNRegion &Region) {
233 return RenderViewOnDC(dc, VPoint);
234}
235
236bool ChartDummy::RenderViewOnDC(wxMemoryDC &dc, const ViewPort &VPoint) {
237 if (m_pBM && m_pBM->IsOk()) {
238 if ((m_pBM->GetWidth() != VPoint.pix_width) ||
239 (m_pBM->GetHeight() != VPoint.pix_height)) {
240 delete m_pBM;
241 m_pBM = NULL;
242 }
243 } else {
244 delete m_pBM;
245 m_pBM = NULL;
246 }
247
248 if (VPoint.pix_width && VPoint.pix_height) {
249 if (NULL == m_pBM)
250 m_pBM = new wxBitmap(VPoint.pix_width, VPoint.pix_height, -1);
251
252 dc.SelectObject(*m_pBM);
253
254 dc.SetBackground(*wxBLACK_BRUSH);
255 dc.Clear();
256 }
257
258 return true;
259}
260
261bool ChartDummy::AdjustVP(ViewPort &vp_last, ViewPort &vp_proposed) {
262 return false;
263}
264
265void ChartDummy::GetValidCanvasRegion(const ViewPort &VPoint,
266 OCPNRegion *pValidRegion) {
267 pValidRegion->Clear();
268 pValidRegion->Union(0, 0, 1, 1);
269}
270
271LLRegion ChartDummy::GetValidRegion() { return LLRegion(); }
272
273// ============================================================================
274// ChartGEO implementation
275// ============================================================================
276ChartGEO::ChartGEO() { m_ChartType = CHART_TYPE_GEO; }
277
278ChartGEO::~ChartGEO() {}
279
280InitReturn ChartGEO::Init(const wxString &name, ChartInitFlag init_flags) {
281#define BUF_LEN_MAX 4096
282
283 PreInit(name, init_flags, GLOBAL_COLOR_SCHEME_DAY);
284
285 char buffer[BUF_LEN_MAX];
286
287 ifs_hdr =
288 new wxFFileInputStream(name); // open the file as a read-only stream
289
290 m_filesize = wxFileName::GetSize(name);
291
292 if (!ifs_hdr->IsOk()) return INIT_FAIL_REMOVE;
293
294 int nPlypoint = 0;
295 Plypoint *pPlyTable = (Plypoint *)malloc(sizeof(Plypoint));
296
297 m_FullPath = name;
298 m_Description = m_FullPath;
299
300 wxFileName GEOFile(m_FullPath);
301
302 wxString Path;
303 Path = GEOFile.GetPath(wxPATH_GET_SEPARATOR | wxPATH_GET_VOLUME);
304
305 // Read the GEO file, extracting useful information
306
307 ifs_hdr->SeekI(0, wxFromStart); // rewind
308
309 Size_X = Size_Y = 0;
310
311 while ((ReadBSBHdrLine(ifs_hdr, &buffer[0], BUF_LEN_MAX)) != 0) {
312 wxString str_buf(buffer, wxConvUTF8);
313 if (!strncmp(buffer, "Bitmap", 6)) {
314 wxStringTokenizer tkz(str_buf, "=");
315 wxString token = tkz.GetNextToken();
316 if (token.IsSameAs("Bitmap", TRUE)) {
317 pBitmapFilePath = new wxString();
318
319 int i;
320 i = tkz.GetPosition();
321 pBitmapFilePath->Clear();
322 while (buffer[i]) {
323 pBitmapFilePath->Append(buffer[i]);
324 i++;
325 }
326 }
327 }
328
329 else if (!strncmp(buffer, "Scale", 5)) {
330 wxStringTokenizer tkz(str_buf, "=");
331 wxString token = tkz.GetNextToken();
332 if (token.IsSameAs("Scale", TRUE)) // extract Scale
333 {
334 int i;
335 i = tkz.GetPosition();
336 m_Chart_Scale = atoi(&buffer[i]);
337 }
338 }
339
340 else if (!strncmp(buffer, "Depth", 5)) {
341 wxStringTokenizer tkz(str_buf, "=");
342 wxString token = tkz.GetNextToken();
343 if (token.IsSameAs("Depth Units", FALSE)) // extract Depth Units
344 {
345 int i;
346 i = tkz.GetPosition();
347 wxString str(&buffer[i], wxConvUTF8);
348 m_DepthUnits = str.Trim();
349 }
350 }
351
352 else if (!strncmp(buffer, "Point", 5)) // Extract RefPoints
353 {
354 int i, xr, yr;
355 float ltr, lnr;
356 sscanf(&buffer[0], "Point%d=%f %f %d %d", &i, &lnr, &ltr, &yr, &xr);
357 pRefTable =
358 (Refpoint *)realloc(pRefTable, sizeof(Refpoint) * (nRefpoint + 1));
359 pRefTable[nRefpoint].xr = xr;
360 pRefTable[nRefpoint].yr = yr;
361 pRefTable[nRefpoint].latr = ltr;
362 pRefTable[nRefpoint].lonr = lnr;
363 pRefTable[nRefpoint].bXValid = 1;
364 pRefTable[nRefpoint].bYValid = 1;
365
366 nRefpoint++;
367
368 }
369
370 else if (!strncmp(buffer, "Vertex", 6)) {
371 int i;
372 float ltp, lnp;
373 sscanf(buffer, "Vertex%d=%f %f", &i, &ltp, &lnp);
374 Plypoint *tmp = pPlyTable;
375 pPlyTable =
376 (Plypoint *)realloc(pPlyTable, sizeof(Plypoint) * (nPlypoint + 1));
377 if (NULL == pPlyTable) {
378 free(tmp);
379 tmp = NULL;
380 } else {
381 pPlyTable[nPlypoint].ltp = ltp;
382 pPlyTable[nPlypoint].lnp = lnp;
383 nPlypoint++;
384 }
385 }
386
387 else if (!strncmp(buffer, "Date Pub", 8)) {
388 char date_string[40];
389 char date_buf[10];
390 sscanf(buffer, "Date Published=%s\r\n", &date_string[0]);
391 wxString date_wxstr(date_string, wxConvUTF8);
392 wxDateTime dt;
393 if (dt.ParseDate(date_wxstr)) // successful parse?
394 {
395 sprintf(date_buf, "%d", dt.GetYear());
396 } else {
397 sscanf(date_string, "%s", date_buf);
398 }
399 m_PubYear = wxString(date_buf, wxConvUTF8);
400 }
401
402 else if (!strncmp(buffer, "Skew", 4)) {
403 wxStringTokenizer tkz(str_buf, "=");
404 wxString token = tkz.GetNextToken();
405 if (token.IsSameAs("Skew Angle", FALSE)) // extract Skew Angle
406 {
407 int i;
408 i = tkz.GetPosition();
409 float fcs;
410 sscanf(&buffer[i], "%f,", &fcs);
411 m_Chart_Skew = fcs;
412 }
413 }
414
415 else if (!strncmp(buffer, "Latitude Offset", 15)) {
416 wxStringTokenizer tkz(str_buf, "=");
417 wxString token = tkz.GetNextToken();
418 if (token.IsSameAs("Latitude Offset", FALSE)) {
419 int i;
420 i = tkz.GetPosition();
421 float lto;
422 sscanf(&buffer[i], "%f,", &lto);
423 m_dtm_lat = lto;
424 }
425 }
426
427 else if (!strncmp(buffer, "Longitude Offset", 16)) {
428 wxStringTokenizer tkz(str_buf, "=");
429 wxString token = tkz.GetNextToken();
430 if (token.IsSameAs("Longitude Offset", FALSE)) {
431 int i;
432 i = tkz.GetPosition();
433 float lno;
434 sscanf(&buffer[i], "%f,", &lno);
435 m_dtm_lon = lno;
436 }
437 }
438
439 else if (!strncmp(buffer, "Datum", 5)) {
440 wxStringTokenizer tkz(str_buf, "=");
441 wxString token = tkz.GetNextToken();
442 if (token.IsSameAs("Datum", FALSE)) {
443 token = tkz.GetNextToken();
444 m_datum_str = token;
445 }
446 }
447
448 else if (!strncmp(buffer, "Name", 4)) {
449 wxStringTokenizer tkz(str_buf, "=");
450 wxString token = tkz.GetNextToken();
451 if (token.IsSameAs("Name", FALSE)) // Name
452 {
453 int i;
454 i = tkz.GetPosition();
455 m_Name.Clear();
456 while (isprint(buffer[i]) && (i < 80)) m_Name.Append(buffer[i++]);
457 }
458 }
459 } // while
460
461 // Extract the remaining data from .NOS Bitmap file
462 ifs_bitmap = NULL;
463
464 // Something wrong with the .geo file, there is no Bitmap reference
465 // This is where the arbitrarily bad file is caught, such as
466 // a file with.GEO extension that is not really a chart
467
468 if (pBitmapFilePath == NULL) {
469 free(pPlyTable);
470 return INIT_FAIL_REMOVE;
471 }
472
473 wxString NOS_Name(*pBitmapFilePath); // take a copy
474
475 wxDir target_dir(Path);
476 wxArrayString file_array;
477 int nfiles = wxDir::GetAllFiles(Path, &file_array);
478 int ifile;
479
480 pBitmapFilePath->Prepend(Path);
481
482 wxFileName NOS_filename(*pBitmapFilePath);
483 if (!NOS_filename.FileExists()) {
484 // File as fetched verbatim from the .geo file doesn't exist.
485 // Try all possible upper/lower cases
486 // Extract the filename and extension
487 wxString fname(NOS_filename.GetName());
488 wxString fext(NOS_filename.GetExt());
489
490 // Try all four combinations, the hard way
491 // case 1
492 fname.MakeLower();
493 fext.MakeLower();
494 NOS_filename.SetName(fname);
495 NOS_filename.SetExt(fext);
496
497 if (NOS_filename.FileExists()) goto found_uclc_file;
498
499 // case 2
500 fname.MakeLower();
501 fext.MakeUpper();
502 NOS_filename.SetName(fname);
503 NOS_filename.SetExt(fext);
504
505 if (NOS_filename.FileExists()) goto found_uclc_file;
506
507 // case 3
508 fname.MakeUpper();
509 fext.MakeLower();
510 NOS_filename.SetName(fname);
511 NOS_filename.SetExt(fext);
512
513 if (NOS_filename.FileExists()) goto found_uclc_file;
514
515 // case 4
516 fname.MakeUpper();
517 fext.MakeUpper();
518 NOS_filename.SetName(fname);
519 NOS_filename.SetExt(fext);
520
521 if (NOS_filename.FileExists()) goto found_uclc_file;
522
523 // Search harder
524
525 for (ifile = 0; ifile < nfiles; ifile++) {
526 wxString file_up = file_array[ifile];
527 file_up.MakeUpper();
528
529 wxString target_up = *pBitmapFilePath;
530 target_up.MakeUpper();
531
532 if (file_up.IsSameAs(target_up)) {
533 NOS_filename.Clear();
534 NOS_filename.Assign(file_array[ifile]);
535 goto found_uclc_file;
536 }
537 }
538
539 free(pPlyTable);
540 return INIT_FAIL_REMOVE; // not found at all
541
542 found_uclc_file:
543
544 delete pBitmapFilePath; // fix up the member element
545 pBitmapFilePath = new wxString(NOS_filename.GetFullPath());
546 }
547 ifss_bitmap =
548 new wxFFileInputStream(*pBitmapFilePath); // open the bitmap file
549 ifs_bitmap = new wxBufferedInputStream(*ifss_bitmap);
550
551 if (!ifss_bitmap->IsOk()) {
552 free(pPlyTable);
553 return INIT_FAIL_REMOVE;
554 }
555
556 while ((ReadBSBHdrLine(ifss_bitmap, &buffer[0], BUF_LEN_MAX)) != 0) {
557 wxString str_buf(buffer, wxConvUTF8);
558
559 if (!strncmp(buffer, "NOS", 3)) {
560 wxStringTokenizer tkz(str_buf, ",=");
561 while (tkz.HasMoreTokens()) {
562 wxString token = tkz.GetNextToken();
563 if (token.IsSameAs("RA", TRUE)) // extract RA=x,y
564 {
565 int i;
566 tkz.GetNextToken();
567 tkz.GetNextToken();
568 i = tkz.GetPosition();
569 Size_X = atoi(&buffer[i]);
570 wxString token = tkz.GetNextToken();
571 i = tkz.GetPosition();
572 Size_Y = atoi(&buffer[i]);
573 } else if (token.IsSameAs("DU", TRUE)) // extract DU=n
574 {
575 token = tkz.GetNextToken();
576 long temp_du;
577 if (token.ToLong(&temp_du)) m_Chart_DU = temp_du;
578 }
579 }
580
581 }
582
583 else if (!strncmp(buffer, "RGB", 3))
584 CreatePaletteEntry(buffer, COLOR_RGB_DEFAULT);
585
586 else if (!strncmp(buffer, "DAY", 3))
587 CreatePaletteEntry(buffer, DAY);
588
589 else if (!strncmp(buffer, "DSK", 3))
590 CreatePaletteEntry(buffer, DUSK);
591
592 else if (!strncmp(buffer, "NGT", 3))
593 CreatePaletteEntry(buffer, NIGHT);
594
595 else if (!strncmp(buffer, "NGR", 3))
596 CreatePaletteEntry(buffer, NIGHTRED);
597
598 else if (!strncmp(buffer, "GRY", 3))
599 CreatePaletteEntry(buffer, GRAY);
600
601 else if (!strncmp(buffer, "PRC", 3))
602 CreatePaletteEntry(buffer, PRC);
603
604 else if (!strncmp(buffer, "PRG", 3))
605 CreatePaletteEntry(buffer, PRG);
606 }
607
608 // Validate some of the header data
609 if (Size_X <= 0 || Size_Y <= 0) {
610 free(pPlyTable);
611 return INIT_FAIL_REMOVE;
612 }
613
614 if (nPlypoint < 3) {
615 wxString msg(" Chart File contains less than 3 PLY points: ");
616 msg.Append(m_FullPath);
617 wxLogMessage(msg);
618 free(pPlyTable);
619
620 return INIT_FAIL_REMOVE;
621 }
622
623 if (m_datum_str.IsEmpty()) {
624 wxString msg(" Chart datum not specified on chart ");
625 msg.Append(m_FullPath);
626 wxLogMessage(msg);
627 wxLogMessage(" Default datum (WGS84) substituted.");
628
629 // return INIT_FAIL_REMOVE;
630 } else {
631 char d_str[100];
632 strncpy(d_str, m_datum_str.mb_str(), 99);
633 d_str[99] = 0;
634
635 int datum_index = GetDatumIndex(d_str);
636
637 if (datum_index < 0) {
638 wxString msg(" Chart datum {");
639 msg += m_datum_str;
640 msg += "} invalid on chart ";
641 msg.Append(m_FullPath);
642 wxLogMessage(msg);
643 wxLogMessage(" Default datum (WGS84) substituted.");
644
645 datum_index = DATUM_INDEX_WGS84;
646 }
647 m_datum_index = datum_index;
648 }
649
650 // Convert captured plypoint information into chart COVR structures
651 m_nCOVREntries = 1;
652 m_pCOVRTablePoints = (int *)malloc(sizeof(int));
653 *m_pCOVRTablePoints = nPlypoint;
654 m_pCOVRTable = (float **)malloc(sizeof(float *));
655 *m_pCOVRTable = (float *)malloc(nPlypoint * 2 * sizeof(float));
656 memcpy(*m_pCOVRTable, pPlyTable, nPlypoint * 2 * sizeof(float));
657
658 free(pPlyTable);
659
660 if (!SetMinMax()) return INIT_FAIL_REMOVE; // have to bail here
661
662 AnalyzeSkew();
663
664 if (init_flags == HEADER_ONLY) return INIT_OK;
665
666 // Advance to the data
667 char c;
668 if ((c = ifs_bitmap->GetC()) != 0x1a) {
669 return INIT_FAIL_REMOVE;
670 }
671 if ((c = ifs_bitmap->GetC()) == 0x0d) {
672 if ((c = ifs_bitmap->GetC()) != 0x0a) {
673 return INIT_FAIL_REMOVE;
674 }
675 if ((c = ifs_bitmap->GetC()) != 0x1a) {
676 return INIT_FAIL_REMOVE;
677 }
678 if ((c = ifs_bitmap->GetC()) != 0x00) {
679 return INIT_FAIL_REMOVE;
680 }
681 }
682
683 else if (c != 0x00) {
684 return INIT_FAIL_REMOVE;
685 }
686
687 // Read the Color table bit size
688 nColorSize = ifs_bitmap->GetC();
689 if (nColorSize == wxEOF || nColorSize <= 0 || nColorSize > 7) {
690 wxString msg(" Invalid nColorSize data, corrupt on chart ");
691 msg.Append(m_FullPath);
692 wxLogMessage(msg);
693 return INIT_FAIL_REMOVE;
694 }
695
696 // Perform common post-init actions in ChartBaseBSB
697 InitReturn pi_ret = PostInit();
698 if (pi_ret != INIT_OK)
699 return pi_ret;
700 else
701 return INIT_OK;
702}
703
704// ============================================================================
705// ChartKAP implementation
706// ============================================================================
707
708ChartKAP::ChartKAP() { m_ChartType = CHART_TYPE_KAP; }
709
710ChartKAP::~ChartKAP() {}
711
712InitReturn ChartKAP::Init(const wxString &name, ChartInitFlag init_flags) {
713#define BUF_LEN_MAX 4096
714
715 ifs_hdr = new ChartDataNonSeekableInputStream(
716 name); // open the Header file as a read-only stream
717
718 if (!ifs_hdr->IsOk()) return INIT_FAIL_REMOVE;
719
720 int nPlypoint = 0;
721 Plypoint *pPlyTable = (Plypoint *)malloc(sizeof(Plypoint));
722
723 PreInit(name, init_flags, GLOBAL_COLOR_SCHEME_DAY);
724
725 char buffer[BUF_LEN_MAX];
726
727 m_FullPath = name;
728 m_Description = m_FullPath;
729
730 // Clear georeferencing coefficients
731 for (int icl = 0; icl < 12; icl++) {
732 wpx[icl] = 0;
733 wpy[icl] = 0;
734 pwx[icl] = 0;
735 pwy[icl] = 0;
736 }
737
738 // Validate the BSB header
739 // by reading some characters into a buffer and looking for BSB\ keyword
740
741 unsigned int TestBlockSize = 1999;
742 ifs_hdr->Read(buffer, TestBlockSize);
743
744 if (ifs_hdr->LastRead() != TestBlockSize) {
745 wxString msg;
746 msg.Printf(" Could not read first %d bytes of header for chart file: ",
747 TestBlockSize);
748 msg.Append(name);
749 wxLogMessage(msg);
750 free(pPlyTable);
751 return INIT_FAIL_REMOVE;
752 }
753
754 unsigned int i;
755 for (i = 0; i < TestBlockSize - 4; i++) {
756 // Test for "BSB/"
757 if (buffer[i + 0] == 'B' && buffer[i + 1] == 'S' && buffer[i + 2] == 'B' &&
758 buffer[i + 3] == '/')
759 break;
760
761 // Test for "NOS/"
762 if (buffer[i + 0] == 'N' && buffer[i + 1] == 'O' && buffer[i + 2] == 'S' &&
763 buffer[i + 3] == '/')
764 break;
765 }
766 if (i == TestBlockSize - 4) {
767 wxString msg(" Chart file has no BSB header, cannot Init.");
768 msg.Append(name);
769 wxLogMessage(msg);
770 free(pPlyTable);
771 return INIT_FAIL_REMOVE;
772 }
773
774 // Read and Parse Chart Header, line by line
775 ifs_hdr->SeekI(0, wxFromStart); // rewind
776
777 Size_X = Size_Y = 0;
778
779 int done_header_parse = 0;
780 wxCSConv iso_conv("ISO-8859-1"); // we will need a converter
781
782 while (done_header_parse == 0) {
783 if (ReadBSBHdrLine(ifs_hdr, buffer, BUF_LEN_MAX) == 0) {
784 unsigned char c;
785 c = ifs_hdr->GetC();
786 ifs_hdr->Ungetch(c);
787
788 if (0x1a == c)
789 done_header_parse = 1;
790 else {
791 free(pPlyTable);
792 return INIT_FAIL_REMOVE;
793 }
794
795 continue;
796 }
797
798 wxString str_buf(buffer, wxConvUTF8);
799 if (!str_buf.Len()) // failed conversion
800 str_buf = wxString(buffer, iso_conv);
801
802 if (str_buf.Find("SHOM") != wxNOT_FOUND) m_b_SHOM = true;
803
804 if (!strncmp(buffer, "BSB", 3)) {
805 wxString clip_str_buf(
806 &buffer[0],
807 iso_conv); // for single byte French encodings of NAme field
808 wxStringTokenizer tkz(clip_str_buf, "/,=");
809 while (tkz.HasMoreTokens()) {
810 wxString token = tkz.GetNextToken();
811 if (token.IsSameAs("RA", TRUE)) // extract RA=x,y
812 {
813 int i;
814 i = tkz.GetPosition();
815 Size_X = atoi(&buffer[i]);
816 wxString token = tkz.GetNextToken();
817 i = tkz.GetPosition();
818 Size_Y = atoi(&buffer[i]);
819 } else if (token.IsSameAs("NA", TRUE)) // extract NA=str
820 {
821 int i = tkz.GetPosition();
822 char nbuf[81];
823 int j = 0;
824 while ((buffer[i] != ',') && (i < 80)) nbuf[j++] = buffer[i++];
825 nbuf[j] = 0;
826 wxString n_str(nbuf, iso_conv);
827 m_Name = n_str;
828 } else if (token.IsSameAs("NU", TRUE)) // extract NU=str
829 {
830 int i = tkz.GetPosition();
831 char nbuf[81];
832 int j = 0;
833 while ((buffer[i] != ',') && (i < 80)) nbuf[j++] = buffer[i++];
834 nbuf[j] = 0;
835 wxString n_str(nbuf, iso_conv);
836 m_ID = n_str;
837 } else if (token.IsSameAs("DU", TRUE)) // extract DU=n
838 {
839 token = tkz.GetNextToken();
840 long temp_du;
841 if (token.ToLong(&temp_du)) m_Chart_DU = temp_du;
842 }
843 }
844 }
845
846 else if (!strncmp(buffer, "KNP", 3)) {
847 wxString conv_buf(buffer, iso_conv);
848 wxStringTokenizer tkz(conv_buf, "/,=");
849 while (tkz.HasMoreTokens()) {
850 wxString token = tkz.GetNextToken();
851 if (token.IsSameAs("SC", TRUE)) // extract Scale
852 {
853 int i;
854 i = tkz.GetPosition();
855 m_Chart_Scale = atoi(&buffer[i]);
856 if (0 == m_Chart_Scale) m_Chart_Scale = 100000000;
857 } else if (token.IsSameAs("SK", TRUE)) // extract Skew
858 {
859 int i;
860 i = tkz.GetPosition();
861 float fcs;
862 sscanf(&buffer[i], "%f,", &fcs);
863 m_Chart_Skew = fcs;
864 } else if (token.IsSameAs("UN", TRUE)) // extract Depth Units
865 {
866 int i;
867 i = tkz.GetPosition();
868 wxString str(&buffer[i], iso_conv);
869 m_DepthUnits = str.BeforeFirst(',');
870 } else if (token.IsSameAs("GD", TRUE)) // extract Datum
871 {
872 int i;
873 i = tkz.GetPosition();
874 wxString str(&buffer[i], iso_conv);
875 m_datum_str = str.BeforeFirst(',').Trim();
876 } else if (token.IsSameAs("SD", TRUE)) // extract Soundings Datum
877 {
878 int i;
879 i = tkz.GetPosition();
880 wxString str(&buffer[i], iso_conv);
881 m_SoundingsDatum = str.BeforeFirst(',').Trim();
882 } else if (token.IsSameAs("PP",
883 TRUE)) // extract Projection Parameter
884 {
885 int i;
886 i = tkz.GetPosition();
887 double fcs;
888 wxString str(&buffer[i], iso_conv);
889 wxString str1 = str.BeforeFirst(',').Trim();
890 if (str1.ToDouble(&fcs)) m_proj_parameter = fcs;
891 } else if (token.IsSameAs("PR", TRUE)) // extract Projection Type
892 {
893 int i;
894 i = tkz.GetPosition();
895 wxString str(&buffer[i], iso_conv);
896 wxString stru = str.MakeUpper();
897 bool bp_set = false;
898 ;
899
900 if (stru.Matches("*MERCATOR*")) {
901 m_projection = PROJECTION_MERCATOR;
902 bp_set = true;
903 }
904
905 if (stru.Matches("*TRANSVERSE*")) {
906 m_projection = PROJECTION_TRANSVERSE_MERCATOR;
907 bp_set = true;
908 }
909
910 if (stru.Matches("*CONIC*")) {
911 m_projection = PROJECTION_POLYCONIC;
912 bp_set = true;
913 }
914
915 if (stru.Matches("*TM*")) {
916 m_projection = PROJECTION_TRANSVERSE_MERCATOR;
917 bp_set = true;
918 }
919
920 if (stru.Matches("*GAUSS CONFORMAL*")) {
921 m_projection = PROJECTION_TRANSVERSE_MERCATOR;
922 bp_set = true;
923 }
924
925 if (!bp_set) {
926 m_projection = PROJECTION_UNKNOWN;
927 wxString msg(" Chart projection is ");
928 msg += tkz.GetNextToken();
929 msg += " which is unsupported. Disabling chart ";
930 msg += m_FullPath;
931 wxLogMessage(msg);
932 free(pPlyTable);
933 return INIT_FAIL_REMOVE;
934 }
935 } else if (token.IsSameAs(
936 "DX",
937 TRUE)) // extract Pixel scale parameter, if present
938 {
939 int i;
940 i = tkz.GetPosition();
941 float x;
942 sscanf(&buffer[i], "%f,", &x);
943 m_dx = x;
944 } else if (token.IsSameAs(
945 "DY",
946 TRUE)) // extract Pixel scale parameter, if present
947 {
948 int i;
949 i = tkz.GetPosition();
950 float x;
951 sscanf(&buffer[i], "%f,", &x);
952 m_dy = x;
953 }
954 }
955 }
956
957 else if (!strncmp(buffer, "RGB", 3))
958 CreatePaletteEntry(buffer, COLOR_RGB_DEFAULT);
959
960 else if (!strncmp(buffer, "DAY", 3))
961 CreatePaletteEntry(buffer, DAY);
962
963 else if (!strncmp(buffer, "DSK", 3))
964 CreatePaletteEntry(buffer, DUSK);
965
966 else if (!strncmp(buffer, "NGT", 3))
967 CreatePaletteEntry(buffer, NIGHT);
968
969 else if (!strncmp(buffer, "NGR", 3))
970 CreatePaletteEntry(buffer, NIGHTRED);
971
972 else if (!strncmp(buffer, "GRY", 3))
973 CreatePaletteEntry(buffer, GRAY);
974
975 else if (!strncmp(buffer, "PRC", 3))
976 CreatePaletteEntry(buffer, PRC);
977
978 else if (!strncmp(buffer, "PRG", 3))
979 CreatePaletteEntry(buffer, PRG);
980
981 else if (!strncmp(buffer, "REF", 3)) {
982 int i, xr, yr;
983 float ltr, lnr;
984 sscanf(&buffer[4], "%d,%d,%d,%f,%f", &i, &xr, &yr, &ltr, &lnr);
985 pRefTable =
986 (Refpoint *)realloc(pRefTable, sizeof(Refpoint) * (nRefpoint + 1));
987 pRefTable[nRefpoint].xr = xr;
988 pRefTable[nRefpoint].yr = yr;
989 pRefTable[nRefpoint].latr = ltr;
990 pRefTable[nRefpoint].lonr = lnr;
991 pRefTable[nRefpoint].bXValid = 1;
992 pRefTable[nRefpoint].bYValid = 1;
993
994 nRefpoint++;
995
996 }
997
998 else if (!strncmp(buffer, "WPX", 3)) {
999 int idx = 0;
1000 double d;
1001 wxStringTokenizer tkz(str_buf.Mid(4), ",");
1002 wxString token = tkz.GetNextToken();
1003
1004 if (token.ToLong((long int *)&wpx_type)) {
1005 while (tkz.HasMoreTokens() && (idx < 12)) {
1006 token = tkz.GetNextToken();
1007 if (token.ToDouble(&d)) {
1008 wpx[idx] = d;
1009 idx++;
1010 }
1011 }
1012 }
1013 n_wpx = idx;
1014 }
1015
1016 else if (!strncmp(buffer, "WPY", 3)) {
1017 int idx = 0;
1018 double d;
1019 wxStringTokenizer tkz(str_buf.Mid(4), ",");
1020 wxString token = tkz.GetNextToken();
1021
1022 if (token.ToLong((long int *)&wpy_type)) {
1023 while (tkz.HasMoreTokens() && (idx < 12)) {
1024 token = tkz.GetNextToken();
1025 if (token.ToDouble(&d)) {
1026 wpy[idx] = d;
1027 idx++;
1028 }
1029 }
1030 }
1031 n_wpy = idx;
1032 }
1033
1034 else if (!strncmp(buffer, "PWX", 3)) {
1035 int idx = 0;
1036 double d;
1037 wxStringTokenizer tkz(str_buf.Mid(4), ",");
1038 wxString token = tkz.GetNextToken();
1039
1040 if (token.ToLong((long int *)&pwx_type)) {
1041 while (tkz.HasMoreTokens() && (idx < 12)) {
1042 token = tkz.GetNextToken();
1043 if (token.ToDouble(&d)) {
1044 pwx[idx] = d;
1045 idx++;
1046 }
1047 }
1048 }
1049 n_pwx = idx;
1050 }
1051
1052 else if (!strncmp(buffer, "PWY", 3)) {
1053 int idx = 0;
1054 double d;
1055 wxStringTokenizer tkz(str_buf.Mid(4), ",");
1056 wxString token = tkz.GetNextToken();
1057
1058 if (token.ToLong((long int *)&pwy_type)) {
1059 while (tkz.HasMoreTokens() && (idx < 12)) {
1060 token = tkz.GetNextToken();
1061 if (token.ToDouble(&d)) {
1062 pwy[idx] = d;
1063 idx++;
1064 }
1065 }
1066 }
1067 n_pwy = idx;
1068 }
1069
1070 else if (!strncmp(buffer, "CPH", 3)) {
1071 float float_cph;
1072 sscanf(&buffer[4], "%f", &float_cph);
1073 m_cph = float_cph;
1074 }
1075
1076 else if (!strncmp(buffer, "VER", 3)) {
1077 wxStringTokenizer tkz(str_buf, "/,=");
1078 wxString token = tkz.GetNextToken();
1079
1080 m_bsb_ver = tkz.GetNextToken();
1081 }
1082
1083 else if (!strncmp(buffer, "DTM", 3)) {
1084 double val;
1085 wxStringTokenizer tkz(str_buf, "/,=");
1086 wxString token = tkz.GetNextToken();
1087
1088 token = tkz.GetNextToken();
1089 if (token.ToDouble(&val)) m_dtm_lat = val;
1090
1091 token = tkz.GetNextToken();
1092 if (token.ToDouble(&val)) m_dtm_lon = val;
1093
1094 // float fdtmlat, fdtmlon;
1095 // sscanf(&buffer[4], "%f,%f", &fdtmlat, &fdtmlon);
1096 // m_dtm_lat = fdtmlat;
1097 // m_dtm_lon = fdtmlon;
1098 }
1099
1100 else if (!strncmp(buffer, "PLY", 3)) {
1101 int i;
1102 float ltp, lnp;
1103 if (sscanf(&buffer[4], "%d,%f,%f", &i, &ltp, &lnp) != 3) {
1104 free(pPlyTable);
1105 return INIT_FAIL_REMOVE;
1106 }
1107 Plypoint *tmp = pPlyTable;
1108 pPlyTable =
1109 (Plypoint *)realloc(pPlyTable, sizeof(Plypoint) * (nPlypoint + 1));
1110 if (NULL == pPlyTable) {
1111 free(tmp);
1112 tmp = NULL;
1113 } else {
1114 pPlyTable[nPlypoint].ltp = ltp;
1115 pPlyTable[nPlypoint].lnp = lnp;
1116 nPlypoint++;
1117 }
1118 if (NULL == pPlyTable || nPlypoint > 1000000) {
1119 // arbitrary 8MB for pPlyTable
1120 nPlypoint = 0;
1121 break;
1122 }
1123 }
1124
1125 else if (!strncmp(buffer, "CED", 3)) {
1126 wxStringTokenizer tkz(str_buf, "/,=");
1127 while (tkz.HasMoreTokens()) {
1128 wxString token = tkz.GetNextToken();
1129 if (token.IsSameAs("ED", TRUE)) // extract Edition Date
1130 {
1131 int i;
1132 i = tkz.GetPosition();
1133
1134 char date_string[40];
1135 char date_buf[16];
1136 date_string[0] = 0;
1137 date_buf[0] = 0;
1138 sscanf(&buffer[i], "%s\r\n", date_string);
1139 wxString date_wxstr(date_string, wxConvUTF8);
1140
1141 wxDateTime dt;
1142 if (dt.ParseDate(date_wxstr)) // successful parse?
1143 {
1144 int iyear =
1145 dt.GetYear(); // GetYear() fails on W98, DMC compiler, wx2.8.3
1146 // BSB charts typically list publish date as xx/yy/zz
1147 // This our own little version of the Y2K problem.
1148 // Just apply some sensible logic
1149
1150 if (iyear < 50) {
1151 iyear += 2000;
1152 dt.SetYear(iyear);
1153 } else if ((iyear >= 50) && (iyear < 100)) {
1154 iyear += 1900;
1155 dt.SetYear(iyear);
1156 }
1157 assert(iyear <= 9999);
1158 sprintf(date_buf, "%d", iyear);
1159
1160 // Initialize the wxDateTime menber for Edition Date
1161 m_EdDate = dt;
1162 } else {
1163 sscanf(date_string, "%s", date_buf);
1164 m_EdDate.Set(1, wxDateTime::Jan,
1165 2000); // Todo this could be smarter
1166 }
1167
1168 m_PubYear = wxString(date_buf, wxConvUTF8);
1169 } else if (token.IsSameAs("SE", TRUE)) // extract Source Edition
1170 {
1171 int i;
1172 i = tkz.GetPosition();
1173 wxString str(&buffer[i], iso_conv);
1174 m_SE = str.BeforeFirst(',');
1175 }
1176 }
1177 }
1178 } // while
1179
1180 // Some charts improperly encode the DTM parameters.
1181 // Identify them as necessary, for further processing
1182 if (m_b_SHOM && (m_bsb_ver == "1.1")) m_b_apply_dtm = false;
1183
1184 // If imbedded coefficients are found,
1185 // then use the polynomial georeferencing algorithms
1186 if (n_pwx && n_pwy && n_pwx && n_pwy) bHaveEmbeddedGeoref = true;
1187
1188 // Set up the projection point according to the projection parameter
1189 if (m_projection == PROJECTION_MERCATOR)
1190 m_proj_lat = m_proj_parameter;
1191 else if (m_projection == PROJECTION_TRANSVERSE_MERCATOR)
1192 m_proj_lon = m_proj_parameter;
1193 else if (m_projection == PROJECTION_POLYCONIC)
1194 m_proj_lon = m_proj_parameter;
1195
1196 // We have seen improperly coded charts, with non-sense value of PP
1197 // parameter FS#1251 Check and override if necessary
1198 if (m_proj_lat > 82.0 || m_proj_lat < -82.0) m_proj_lat = 0.0;
1199
1200 // Validate some of the header data
1201 if (Size_X <= 0 || Size_Y <= 0) {
1202 free(pPlyTable);
1203 return INIT_FAIL_REMOVE;
1204 }
1205
1206 if (nPlypoint < 3) {
1207 wxString msg(" Chart File contains less than 3 or too many PLY points: ");
1208 msg.Append(m_FullPath);
1209 wxLogMessage(msg);
1210 free(pPlyTable);
1211 return INIT_FAIL_REMOVE;
1212 }
1213
1214 if (m_datum_str.IsEmpty()) {
1215 wxString msg(" Chart datum not specified on chart ");
1216 msg.Append(m_FullPath);
1217 wxLogMessage(msg);
1218 wxLogMessage(" Default datum (WGS84) substituted.");
1219
1220 // return INIT_FAIL_REMOVE;
1221 } else {
1222 char d_str[100];
1223 strncpy(d_str, m_datum_str.mb_str(), 99);
1224 d_str[99] = 0;
1225
1226 int datum_index = GetDatumIndex(d_str);
1227
1228 if (datum_index < 0) {
1229 wxString msg(" Chart datum {");
1230 msg += m_datum_str;
1231 msg += "} invalid on chart ";
1232 msg.Append(m_FullPath);
1233 wxLogMessage(msg);
1234 wxLogMessage(" Default datum (WGS84) substituted.");
1235
1236 // return INIT_FAIL_REMOVE;
1237 }
1238 }
1239
1240 /* Augment ply points
1241 This is needed for example on polyconic charts or skewed charts because
1242 straight lines in the chart coordinates can not use simple
1243 interpolation in lat/lon or mercator coordinate space to draw the
1244 borders or be used for quilting operation.
1245 TODO: should this be added as a subroutine for GEO chartso? */
1246 if ((m_projection != PROJECTION_MERCATOR &&
1247 m_projection != PROJECTION_TRANSVERSE_MERCATOR) ||
1248 m_Chart_Skew > 2) {
1249 // Analyze Refpoints early because we need georef coefficient here.
1250 AnalyzeRefpoints(false); // no post test needed
1251
1252 // We need to compute a tentative min/max lat/lon to perform georefs
1253 // These lat/lon extents will be more accurately updated later.
1254 m_LonMax = -360.0;
1255 m_LonMin = 360.0;
1256 m_LatMax = -90.0;
1257 m_LatMin = 90.0;
1258
1259 for (int i = 0; i < nPlypoint; i++) {
1260 m_LatMax = wxMax(m_LatMax, pPlyTable[i].ltp);
1261 m_LatMin = wxMin(m_LatMin, pPlyTable[i].ltp);
1262 m_LonMax = wxMax(m_LonMax, pPlyTable[i].lnp);
1263 m_LonMin = wxMin(m_LonMin, pPlyTable[i].lnp);
1264 }
1265
1266 int count = nPlypoint;
1267 nPlypoint = 0;
1268 Plypoint *pOldPlyTable = pPlyTable;
1269 pPlyTable = NULL;
1270 double lastplylat = 0.0, lastplylon = 0.0, x1 = 0.0, y1 = 0.0, x2, y2;
1271 double plylat, plylon;
1272 for (int i = 0; i < count + 1; i++) {
1273 plylat = pOldPlyTable[i % count].ltp;
1274 plylon = pOldPlyTable[i % count].lnp;
1275 latlong_to_chartpix(plylat, plylon, x2, y2);
1276 if (i > 0) {
1277 if (lastplylon - plylon > 180.)
1278 lastplylon -= 360.;
1279 else if (lastplylon - plylon < -180.)
1280 lastplylon += 360.;
1281
1282 // use 2 degree steps
1283 double steps =
1284 ceil((fabs(lastplylat - plylat) + fabs(lastplylon - plylon)) / 2);
1285 for (double c = 0; c < steps; c++) {
1286 double d = c / steps, lat, lon;
1287 wxPoint2DDouble s;
1288 double x = (1 - d) * x1 + d * x2, y = (1 - d) * y1 + d * y2;
1289 chartpix_to_latlong(x, y, &lat, &lon);
1290 pPlyTable = (Plypoint *)realloc(pPlyTable,
1291 sizeof(Plypoint) * (nPlypoint + 1));
1292 pPlyTable[nPlypoint].ltp = lat;
1293 pPlyTable[nPlypoint].lnp = lon;
1294 nPlypoint++;
1295 }
1296 }
1297 x1 = x2, y1 = y2;
1298 lastplylat = plylat, lastplylon = plylon;
1299 }
1300 free(pOldPlyTable);
1301 }
1302
1303 // Convert captured plypoint information into chart COVR structures
1304
1305 // A special-case test for poorly formatted charts
1306 // We look for cases where the declared PlyPoints are far outside of the
1307 // chart raster bitmap If found, we change the COVR region to the valid
1308 // bitmap region, instead of the default PlyPoints region
1309 // Set a tentative lat/lon range.
1310 m_LonMax = -360.;
1311 m_LonMin = 360.;
1312 for (int i = 0; i < nPlypoint; i++) {
1313 m_LonMin = wxMin(m_LonMin, pPlyTable[i].lnp);
1314 m_LonMax = wxMax(m_LonMax, pPlyTable[i].lnp);
1315 }
1316 // This test does not really work for charts that cross IDL
1317 bool b_test = true;
1318 bool b_adjusted = false;
1319 if (m_LonMax * m_LonMin < 0) {
1320 if ((m_LonMax - m_LonMin) > 180.) b_test = false;
1321 }
1322
1323 if (b_test) {
1324 if (!bHaveEmbeddedGeoref) {
1325 // Analyze Refpoints early because we might need georef coefficient
1326 // here.
1327 AnalyzeRefpoints(false); // no post test needed
1328 }
1329
1330 bool bAdjustPly = false;
1331 wxRect bitRect(0, 0, Size_X, Size_Y);
1332 bitRect.Inflate(5); // allow for a little roundoff error
1333 for (int i = 0; i < nPlypoint; i++) {
1334 double pix_x, pix_y;
1335 latlong_to_chartpix(pPlyTable[i].ltp, pPlyTable[i].lnp, pix_x, pix_y);
1336 if (!bitRect.Contains(pix_x, pix_y)) {
1337 bAdjustPly = true;
1338 if (m_b_cdebug)
1339 printf("Adjusting COVR region on: %s\n", name.ToUTF8().data());
1340 break;
1341 }
1342 }
1343
1344 if (bAdjustPly) {
1345 float *points = new float[2 * nPlypoint];
1346 for (int i = 0; i < nPlypoint; i++)
1347 points[2 * i + 0] = pPlyTable[i].ltp,
1348 points[2 * i + 1] = pPlyTable[i].lnp;
1349 LLRegion covrRegion(nPlypoint, points);
1350 delete[] points;
1351 covrRegion.Intersect(GetValidRegion());
1352
1353 if (covrRegion.contours.size()) { // Check for no intersection caused by
1354 // bogus georef....
1355 m_nCOVREntries = covrRegion.contours.size();
1356 m_pCOVRTablePoints = (int *)malloc(m_nCOVREntries * sizeof(int));
1357 m_pCOVRTable = (float **)malloc(m_nCOVREntries * sizeof(float *));
1358 std::list<poly_contour>::iterator it = covrRegion.contours.begin();
1359 for (int i = 0; i < m_nCOVREntries; i++) {
1360 m_pCOVRTablePoints[i] = it->size();
1361 m_pCOVRTable[i] =
1362 (float *)malloc(m_pCOVRTablePoints[i] * 2 * sizeof(float));
1363 std::list<contour_pt>::iterator jt = it->begin();
1364 for (int j = 0; j < m_pCOVRTablePoints[i]; j++) {
1365 m_pCOVRTable[i][2 * j + 0] = jt->y;
1366 m_pCOVRTable[i][2 * j + 1] = jt->x;
1367 jt++;
1368 }
1369 it++;
1370 }
1371 b_adjusted = true;
1372 }
1373 }
1374 }
1375
1376 if (!b_adjusted) {
1377 m_nCOVREntries = 1;
1378 m_pCOVRTablePoints = (int *)malloc(sizeof(int));
1379 *m_pCOVRTablePoints = nPlypoint;
1380 m_pCOVRTable = (float **)malloc(sizeof(float *));
1381 *m_pCOVRTable = (float *)malloc(nPlypoint * 2 * sizeof(float));
1382 memcpy(*m_pCOVRTable, pPlyTable, nPlypoint * 2 * sizeof(float));
1383 }
1384
1385 free(pPlyTable);
1386
1387 // Setup the datum transform parameters
1388 char d_str[100];
1389 strncpy(d_str, m_datum_str.mb_str(), 99);
1390 d_str[99] = 0;
1391
1392 int datum_index = GetDatumIndex(d_str);
1393 m_datum_index = datum_index;
1394
1395 if (datum_index < 0)
1396 m_ExtraInfo = "---<<< Warning: Chart Datum may be incorrect. >>>---";
1397
1398 // Establish defaults, may be overridden later
1399 m_lon_datum_adjust = (-m_dtm_lon) / 3600.;
1400 m_lat_datum_adjust = (-m_dtm_lat) / 3600.;
1401
1402 // Adjust the PLY points to WGS84 datum
1403 Plypoint *ppp = (Plypoint *)GetCOVRTableHead(0);
1404 int cnPlypoint = GetCOVRTablenPoints(0);
1405
1406 for (int u = 0; u < cnPlypoint; u++) {
1407 double dlat = 0;
1408 double dlon = 0;
1409
1410 if (m_datum_index == DATUM_INDEX_WGS84 ||
1411 m_datum_index == DATUM_INDEX_UNKNOWN) {
1412 dlon = m_dtm_lon / 3600.;
1413 dlat = m_dtm_lat / 3600.;
1414 }
1415
1416 else {
1417 double to_lat, to_lon;
1418 MolodenskyTransform(ppp->ltp, ppp->lnp, &to_lat, &to_lon, m_datum_index,
1419 DATUM_INDEX_WGS84);
1420 dlon = (to_lon - ppp->lnp);
1421 dlat = (to_lat - ppp->ltp);
1422 if (m_b_apply_dtm) {
1423 dlon += m_dtm_lon / 3600.;
1424 dlat += m_dtm_lat / 3600.;
1425 }
1426 }
1427
1428 ppp->lnp += dlon;
1429 ppp->ltp += dlat;
1430 ppp++;
1431 }
1432
1433 if (!SetMinMax()) return INIT_FAIL_REMOVE; // have to bail here
1434
1435 AnalyzeSkew();
1436
1437 if (init_flags == HEADER_ONLY) return INIT_OK;
1438
1439 // Advance to the data
1440 unsigned char c;
1441 bool bcorrupt = false;
1442
1443 if ((c = ifs_hdr->GetC()) != 0x1a) {
1444 bcorrupt = true;
1445 }
1446 if ((c = ifs_hdr->GetC()) == 0x0d) {
1447 if ((c = ifs_hdr->GetC()) != 0x0a) {
1448 bcorrupt = true;
1449 }
1450 if ((c = ifs_hdr->GetC()) != 0x1a) {
1451 bcorrupt = true;
1452 }
1453 if ((c = ifs_hdr->GetC()) != 0x00) {
1454 bcorrupt = true;
1455 }
1456 }
1457
1458 else if (c != 0x00) {
1459 bcorrupt = true;
1460 }
1461
1462 if (bcorrupt) {
1463 wxString msg(" Chart File RLL data corrupt on chart ");
1464 msg.Append(m_FullPath);
1465 wxLogMessage(msg);
1466
1467 return INIT_FAIL_REMOVE;
1468 }
1469
1470 // Read the Color table bit size
1471 nColorSize = ifs_hdr->GetC();
1472 if (nColorSize == wxEOF || nColorSize <= 0 || nColorSize > 7) {
1473 wxString msg(" Invalid nColorSize data, corrupt on chart ");
1474 msg.Append(m_FullPath);
1475 wxLogMessage(msg);
1476 return INIT_FAIL_REMOVE;
1477 }
1478
1479 nFileOffsetDataStart = ifs_hdr->TellI();
1480 delete ifs_hdr;
1481 ifs_hdr = NULL;
1482
1483 ChartDataInputStream *stream =
1484 new ChartDataInputStream(name); // Open again, as the bitmap
1485 wxString tempfile;
1486#ifdef OCPN_USE_LZMA
1487 tempfile = stream->TempFileName();
1488#endif
1489 m_filesize = wxFileName::GetSize(tempfile.empty() ? name : tempfile);
1490
1491 ifss_bitmap = stream;
1492 ifs_bitmap = new wxBufferedInputStream(*ifss_bitmap);
1493
1494 // Perform common post-init actions in ChartBaseBSB
1495 InitReturn pi_ret = PostInit();
1496 if (pi_ret != INIT_OK) return pi_ret;
1497 return INIT_OK;
1498}
1499
1500// ============================================================================
1501// ChartBaseBSB implementation
1502// ============================================================================
1503
1504ChartBaseBSB::ChartBaseBSB() {
1505 // Init some private data
1506 m_ChartFamily = CHART_FAMILY_RASTER;
1507
1508 pBitmapFilePath = NULL;
1509
1510 pline_table = NULL;
1511 ifs_buf = NULL;
1512
1513 cached_image_ok = 0;
1514
1515 pRefTable = (Refpoint *)malloc(sizeof(Refpoint));
1516 nRefpoint = 0;
1517 cPoints.status = 0;
1518 bHaveEmbeddedGeoref = false;
1519 n_wpx = 0;
1520 n_wpy = 0;
1521 n_pwx = 0;
1522 n_pwy = 0;
1523
1524#ifdef __ANDROID__
1525 bUseLineCache = false;
1526#else
1527 bUseLineCache = true;
1528#endif
1529
1530 m_Chart_Skew = 0.0;
1531
1532 pPixCache = NULL;
1533
1534 pLineCache = NULL;
1535
1536 m_bilinear_limit = 8; // bilinear scaling only up to n
1537
1538 ifs_bitmap = NULL;
1539 ifss_bitmap = NULL;
1540 ifs_hdr = NULL;
1541
1542 for (int i = 0; i < N_BSB_COLORS; i++) pPalettes[i] = NULL;
1543
1544 bGeoErrorSent = false;
1545 m_Chart_DU = 0;
1546 m_cph = 0.;
1547
1548 m_mapped_color_index = COLOR_RGB_DEFAULT;
1549
1550 m_datum_str = "WGS84"; // assume until proven otherwise
1551
1552 m_dtm_lat = 0.;
1553 m_dtm_lon = 0.;
1554
1555 m_dx = 0.;
1556 m_dy = 0.;
1557 m_proj_lat = 0.;
1558 m_proj_lon = 0.;
1559 m_proj_parameter = 0.;
1560 m_b_SHOM = false;
1561 m_b_apply_dtm = true;
1562
1563 m_b_cdebug = 0;
1564
1565#ifdef OCPN_USE_CONFIG
1566 wxFileConfig *pfc = (wxFileConfig *)pConfig;
1567 pfc->SetPath("/Settings");
1568 pfc->Read("DebugBSBImg", &m_b_cdebug, 0);
1569#endif
1570}
1571
1572ChartBaseBSB::~ChartBaseBSB() {
1573 if (pBitmapFilePath) delete pBitmapFilePath;
1574
1575 if (pline_table) free(pline_table);
1576
1577 if (ifs_buf) free(ifs_buf);
1578
1579 free(pRefTable);
1580 // free(pPlyTable);
1581
1582 delete ifs_bitmap;
1583 delete ifs_hdr;
1584 delete ifss_bitmap;
1585
1586 if (cPoints.status) {
1587 free(cPoints.tx);
1588 free(cPoints.ty);
1589 free(cPoints.lon);
1590 free(cPoints.lat);
1591
1592 free(cPoints.pwx);
1593 free(cPoints.wpx);
1594 free(cPoints.pwy);
1595 free(cPoints.wpy);
1596 }
1597
1598 // Free the line cache
1599 FreeLineCacheRows();
1600 free(pLineCache);
1601
1602 delete pPixCache;
1603
1604 for (int i = 0; i < N_BSB_COLORS; i++) delete pPalettes[i];
1605}
1606
1607void ChartBaseBSB::FreeLineCacheRows(int start, int end) {
1608 if (pLineCache) {
1609 if (end < 0)
1610 end = Size_Y;
1611 else
1612 end = wxMin(end, Size_Y);
1613 for (int ylc = start; ylc < end; ylc++) {
1614 CachedLine *pt = &pLineCache[ylc];
1615 if (pt->bValid) {
1616 free(pt->pTileOffset);
1617 free(pt->pPix);
1618 pt->bValid = false;
1619 }
1620 }
1621 }
1622}
1623
1624bool ChartBaseBSB::HaveLineCacheRow(int row) {
1625 if (pLineCache) {
1626 CachedLine *pt = &pLineCache[row];
1627 return pt->bValid;
1628 }
1629 return false;
1630}
1631
1632// Report recommended minimum and maximum scale values for which use of this
1633// chart is valid
1634
1635double ChartBaseBSB::GetNormalScaleMin(double canvas_scale_factor,
1636 bool b_allow_overzoom) {
1637 // if(b_allow_overzoom)
1638 return (canvas_scale_factor / m_ppm_avg) /
1639 32; // allow wide range overzoom overscale
1640 // else
1641 // return (canvas_scale_factor / m_ppm_avg) / 2; // don't
1642 // suggest too much overscale
1643}
1644
1645double ChartBaseBSB::GetNormalScaleMax(double canvas_scale_factor,
1646 int canvas_width) {
1647 return (canvas_scale_factor / m_ppm_avg) *
1648 4.0; // excessive underscale is slow, and unreadable
1649}
1650
1651double ChartBaseBSB::GetNearestPreferredScalePPM(double target_scale_ppm) {
1653 target_scale_ppm, .01, 64.); // changed from 32 to 64 to allow super
1654 // small scale BSB charts as quilt base
1655}
1656
1658 double scale_factor_min,
1659 double scale_factor_max) {
1660 double chart_1x_scale = GetPPM();
1661
1662 double binary_scale_factor = 1.;
1663
1664 // Overzoom....
1665 if (chart_1x_scale > target_scale) {
1666 double binary_scale_factor_max = 1 / scale_factor_min;
1667
1668 while (binary_scale_factor < binary_scale_factor_max) {
1669 if (fabs((chart_1x_scale / binary_scale_factor) - target_scale) <
1670 (target_scale * 0.05))
1671 break;
1672 if ((chart_1x_scale / binary_scale_factor) < target_scale)
1673 break;
1674 else
1675 binary_scale_factor *= 2.;
1676 }
1677 }
1678
1679 // Underzoom.....
1680 else {
1681 int ibsf = 1;
1682 int isf_max = (int)scale_factor_max;
1683 while (ibsf < isf_max) {
1684 if (fabs((chart_1x_scale * ibsf) - target_scale) < (target_scale * 0.05))
1685 break;
1686
1687 else if ((chart_1x_scale * ibsf) > target_scale) {
1688 if (ibsf > 1) ibsf /= 2;
1689 break;
1690 } else
1691 ibsf *= 2;
1692 }
1693
1694 binary_scale_factor = 1. / ibsf;
1695 }
1696
1697 return chart_1x_scale / binary_scale_factor;
1698}
1699
1700InitReturn ChartBaseBSB::Init(const wxString &name, ChartInitFlag init_flags) {
1701 m_global_color_scheme = GLOBAL_COLOR_SCHEME_RGB;
1702 return INIT_OK;
1703}
1704
1705InitReturn ChartBaseBSB::PreInit(const wxString &name, ChartInitFlag init_flags,
1706 ColorScheme cs) {
1707 m_global_color_scheme = cs;
1708 return INIT_OK;
1709}
1710
1711void ChartBaseBSB::CreatePaletteEntry(char *buffer, int palette_index) {
1712 if (palette_index < N_BSB_COLORS) {
1713 if (!pPalettes[palette_index]) pPalettes[palette_index] = new opncpnPalette;
1714 opncpnPalette *pp = pPalettes[palette_index];
1715
1716 pp->FwdPalette =
1717 (int *)realloc(pp->FwdPalette, (pp->nFwd + 1) * sizeof(int));
1718 pp->RevPalette =
1719 (int *)realloc(pp->RevPalette, (pp->nRev + 1) * sizeof(int));
1720 pp->nFwd++;
1721 pp->nRev++;
1722
1723 int i;
1724 int n, r, g, b;
1725 sscanf(&buffer[4], "%d,%d,%d,%d", &n, &r, &g, &b);
1726
1727 i = n;
1728
1729 int fcolor, rcolor;
1730 fcolor = (b << 16) + (g << 8) + r;
1731 rcolor = (r << 16) + (g << 8) + b;
1732
1733 pp->RevPalette[i] = rcolor;
1734 pp->FwdPalette[i] = fcolor;
1735 }
1736}
1737
1738InitReturn ChartBaseBSB::PostInit() {
1739 // catch undefined shift if not already done in derived classes
1740 if (nColorSize == wxEOF || nColorSize <= 0 || nColorSize > 7) {
1741 wxString msg(" Invalid nColorSize data, corrupt in PostInit() on chart ");
1742 msg.Append(m_FullPath);
1743 wxLogMessage(msg);
1744 return INIT_FAIL_REMOVE;
1745 }
1746
1747 if (Size_X <= 0 || Size_X > INT_MAX / 4 || Size_Y <= 0 ||
1748 Size_Y - 1 > INT_MAX / 4) {
1749 wxString msg(
1750 " Invalid Size_X/Size_Y data, corrupt in PostInit() on chart ");
1751 msg.Append(m_FullPath);
1752 wxLogMessage(msg);
1753 return INIT_FAIL_REMOVE;
1754 }
1755
1756 // Validate the palette array, substituting DEFAULT for missing entries
1757 int nfwd_def = 1;
1758 int nrev_def = 1;
1759 if (pPalettes[COLOR_RGB_DEFAULT]) {
1760 nrev_def = pPalettes[COLOR_RGB_DEFAULT]->nRev;
1761 nfwd_def = pPalettes[COLOR_RGB_DEFAULT]->nFwd;
1762 }
1763
1764 for (int i = 0; i < N_BSB_COLORS; i++) {
1765 if (pPalettes[i] == NULL) {
1766 opncpnPalette *pNullSubPal = new opncpnPalette;
1767
1768 pNullSubPal->nFwd = nfwd_def; // copy the palette count
1769 pNullSubPal->nRev = nrev_def; // copy the palette count
1770 // Deep copy the palette rgb tables
1771 free(pNullSubPal->FwdPalette);
1772 pNullSubPal->FwdPalette = (int *)malloc(pNullSubPal->nFwd * sizeof(int));
1773 if (pPalettes[COLOR_RGB_DEFAULT])
1774 memcpy(pNullSubPal->FwdPalette,
1775 pPalettes[COLOR_RGB_DEFAULT]->FwdPalette,
1776 pNullSubPal->nFwd * sizeof(int));
1777
1778 free(pNullSubPal->RevPalette);
1779 pNullSubPal->RevPalette = (int *)malloc(pNullSubPal->nRev * sizeof(int));
1780 if (pPalettes[COLOR_RGB_DEFAULT])
1781 memcpy(pNullSubPal->RevPalette,
1782 pPalettes[COLOR_RGB_DEFAULT]->RevPalette,
1783 pNullSubPal->nRev * sizeof(int));
1784
1785 pPalettes[i] = pNullSubPal;
1786 }
1787 }
1788
1789 // Establish the palette type and default palette
1790 palette_direction = GetPaletteDir();
1791
1792 SetColorScheme(m_global_color_scheme, false);
1793
1794 // Allocate memory for ifs file buffering
1795 ifs_bufsize = Size_X * 4;
1796 ifs_buf = (unsigned char *)malloc(ifs_bufsize);
1797 if (!ifs_buf) return INIT_FAIL_REMOVE;
1798
1799 ifs_bufend = ifs_buf + ifs_bufsize;
1800 ifs_lp = ifs_bufend;
1801 ifs_file_offset = -ifs_bufsize;
1802
1803 // Create and load the line offset index table
1804 pline_table = NULL;
1805 pline_table = (int *)malloc((Size_Y + 1) * sizeof(int)); // Ugly....
1806 if (!pline_table) return INIT_FAIL_REMOVE;
1807
1808 ifs_bitmap->SeekI((Size_Y + 1) * -4,
1809 wxFromEnd); // go to Beginning of offset table
1810 pline_table[Size_Y] = ifs_bitmap->TellI(); // fill in useful last table entry
1811
1812 unsigned char *tmp = (unsigned char *)malloc(Size_Y * sizeof(int));
1813 ifs_bitmap->Read(tmp, Size_Y * sizeof(int));
1814 if (ifs_bitmap->LastRead() != Size_Y * sizeof(int)) {
1815 wxString msg(" Chart File corrupt in PostInit() on chart ");
1816 msg.Append(m_FullPath);
1817 wxLogMessage(msg);
1818 free(tmp);
1819
1820 return INIT_FAIL_REMOVE;
1821 }
1822
1823 int offset;
1824 unsigned char *b = tmp;
1825 for (int ifplt = 0; ifplt < Size_Y; ifplt++) {
1826 offset = 0;
1827 offset += *b++ * 256 * 256 * 256;
1828 offset += *b++ * 256 * 256;
1829 offset += *b++ * 256;
1830 offset += *b++;
1831
1832 pline_table[ifplt] = offset;
1833 }
1834 free(tmp);
1835 // Try to validate the line index
1836
1837 bool bline_index_ok = true;
1838 m_nLineOffset = 0;
1839
1840 wxULongLong bitmap_filesize = m_filesize;
1841 if ((m_ChartType == CHART_TYPE_GEO) && pBitmapFilePath)
1842 bitmap_filesize = wxFileName::GetSize(*pBitmapFilePath);
1843
1844 // look logically at the line offset table
1845 for (int iplt = 0; iplt < Size_Y - 1; iplt++) {
1846 if (pline_table[iplt] > bitmap_filesize) {
1847 wxString msg(" Chart File corrupt in PostInit() on chart ");
1848 msg.Append(m_FullPath);
1849 wxLogMessage(msg);
1850
1851 return INIT_FAIL_REMOVE;
1852 }
1853
1854 int thisline_size = pline_table[iplt + 1] - pline_table[iplt];
1855 if (thisline_size < 0) {
1856 wxString msg(" Chart File corrupt in PostInit() on chart ");
1857 msg.Append(m_FullPath);
1858 wxLogMessage(msg);
1859
1860 return INIT_FAIL_REMOVE;
1861 }
1862 }
1863
1864 // For older charts, say Version 1.x, we will try to read the chart and check
1865 // the lines for coherence These older charts are more likely to have index
1866 // troubles.... We only need to check a few lines. Errors are quickly
1867 // apparent.
1868 double ver;
1869 m_bsb_ver.ToDouble(&ver);
1870 if (ver < 2.0) {
1871 for (int iplt = 0; iplt < 10; iplt++) {
1872 if (wxInvalidOffset ==
1873 ifs_bitmap->SeekI(pline_table[iplt], wxFromStart)) {
1874 wxString msg(" Chart File corrupt in PostInit() on chart ");
1875 msg.Append(m_FullPath);
1876 wxLogMessage(msg);
1877
1878 return INIT_FAIL_REMOVE;
1879 }
1880
1881 int thisline_size = pline_table[iplt + 1] - pline_table[iplt];
1882 ifs_bitmap->Read(ifs_buf, thisline_size);
1883
1884 unsigned char *lp = ifs_buf;
1885
1886 unsigned char byNext;
1887 int nLineMarker = 0;
1888 do {
1889 byNext = *lp++;
1890 nLineMarker = nLineMarker * 128 + (byNext & 0x7f);
1891 } while ((byNext & 0x80) != 0);
1892
1893 // Linemarker Correction factor needed here
1894 // Some charts start with LineMarker = 0, some with LineMarker = 1
1895 // Assume the first LineMarker found is the index base, and use
1896 // as a correction offset
1897
1898 if (iplt == 0) m_nLineOffset = nLineMarker;
1899
1900 if (nLineMarker != iplt + m_nLineOffset) {
1901 bline_index_ok = false;
1902 break;
1903 }
1904 }
1905 }
1906
1907 // Recreate the scan line index if the embedded version seems corrupt
1908 if (!bline_index_ok) {
1909 wxString msg(" Line Index corrupt, recreating Index for chart ");
1910 msg.Append(m_FullPath);
1911 wxLogMessage(msg);
1912 if (!CreateLineIndex()) {
1913 wxString msg(" Error creating Line Index for chart ");
1914 msg.Append(m_FullPath);
1915 wxLogMessage(msg);
1916 return INIT_FAIL_REMOVE;
1917 }
1918 }
1919
1920 // Allocate the Line Cache
1921 if (bUseLineCache) {
1922 pLineCache = (CachedLine *)malloc(Size_Y * sizeof(CachedLine));
1923 CachedLine *pt;
1924
1925 for (int ylc = 0; ylc < Size_Y; ylc++) {
1926 pt = &pLineCache[ylc];
1927 pt->bValid = false;
1928 pt->pPix = NULL; //(unsigned char *)malloc(1);
1929 pt->pTileOffset = NULL;
1930 }
1931 } else
1932 pLineCache = NULL;
1933
1934 // Validate/Set Depth Unit Type
1935 wxString test_str = m_DepthUnits.Upper();
1936 if (test_str.IsSameAs("FEET", FALSE))
1937 m_depth_unit_id = DEPTH_UNIT_FEET;
1938 else if (test_str.IsSameAs("METERS", FALSE))
1939 m_depth_unit_id = DEPTH_UNIT_METERS;
1940 else if (test_str.IsSameAs("METRES",
1941 FALSE)) // Special case for alternate spelling
1942 m_depth_unit_id = DEPTH_UNIT_METERS;
1943 else if (test_str.IsSameAs("METRIC", FALSE))
1944 m_depth_unit_id = DEPTH_UNIT_METERS;
1945 else if (test_str.IsSameAs("FATHOMS", FALSE))
1946 m_depth_unit_id = DEPTH_UNIT_FATHOMS;
1947 else if (test_str.Find("FATHOMS") !=
1948 wxNOT_FOUND) // Special case for "Fathoms and Feet"
1949 m_depth_unit_id = DEPTH_UNIT_FATHOMS;
1950 else if (test_str.Find("METERS") !=
1951 wxNOT_FOUND) // Special case for "Meters and decimeters"
1952 m_depth_unit_id = DEPTH_UNIT_METERS;
1953
1954 // Analyze Refpoints
1955 int analyze_ret_val = AnalyzeRefpoints();
1956 if (0 != analyze_ret_val) return INIT_FAIL_REMOVE;
1957
1958 bReadyToRender = true;
1959 return INIT_OK;
1960}
1961
1962bool ChartBaseBSB::CreateLineIndex() {
1963 // Assumes file stream ifs_bitmap is currently open
1964
1965 // wxBufferedInputStream *pbis = new wxBufferedInputStream(*ifss_bitmap);
1966
1967 // Seek to start of data
1968 ifs_bitmap->SeekI(nFileOffsetDataStart); // go to Beginning of data
1969
1970 for (int iplt = 0; iplt < Size_Y; iplt++) {
1971 int offset = ifs_bitmap->TellI();
1972
1973 int iscan = BSBScanScanline(ifs_bitmap);
1974
1975 // There is no sense reporting an error here, since we are recreating after
1976 // an error
1977 /*
1978 if(iscan > Size_Y)
1979 {
1980
1981 wxString msg("CreateLineIndex() failed on chart ");
1982 msg.Append(m_FullPath);
1983 wxLogMessage(msg);
1984 return false;
1985 }
1986
1987 // Skipped lines?
1988 if(iscan != iplt)
1989 {
1990 while((iplt < iscan) && (iplt < Size_Y))
1991 {
1992 pline_table[iplt] = 0;
1993 iplt++;
1994 }
1995 }
1996 */
1997 pline_table[iplt] = offset;
1998 }
1999
2000 return true;
2001}
2002
2003// Invalidate and Free the line cache contents
2004void ChartBaseBSB::InvalidateLineCache() {
2005 if (pLineCache) {
2006 CachedLine *pt;
2007 for (int ylc = 0; ylc < Size_Y; ylc++) {
2008 pt = &pLineCache[ylc];
2009 if (pt) {
2010 free(pt->pPix);
2011 pt->pPix = NULL;
2012 free(pt->pTileOffset);
2013 pt->pTileOffset = NULL;
2014 pt->bValid = false;
2015 }
2016 }
2017 }
2018}
2019
2020bool ChartBaseBSB::GetChartExtent(Extent *pext) {
2021 pext->NLAT = m_LatMax;
2022 pext->SLAT = m_LatMin;
2023 pext->ELON = m_LonMax;
2024 pext->WLON = m_LonMin;
2025
2026 return true;
2027}
2028
2029bool ChartBaseBSB::SetMinMax() {
2030 // Calculate the Chart Extents(M_LatMin, M_LonMin, etc.)
2031 // from the COVR data, for fast database search
2032 m_LonMax = -360.0;
2033 m_LonMin = 360.0;
2034 m_LatMax = -90.0;
2035 m_LatMin = 90.0;
2036
2037 Plypoint *ppp = (Plypoint *)GetCOVRTableHead(0);
2038 int cnPlypoint = GetCOVRTablenPoints(0);
2039
2040 for (int u = 0; u < cnPlypoint; u++) {
2041 if (ppp->lnp > m_LonMax) m_LonMax = ppp->lnp;
2042 if (ppp->lnp < m_LonMin) m_LonMin = ppp->lnp;
2043
2044 if (ppp->ltp > m_LatMax) m_LatMax = ppp->ltp;
2045 if (ppp->ltp < m_LatMin) m_LatMin = ppp->ltp;
2046
2047 ppp++;
2048 }
2049
2050 // Check for special cases
2051
2052 // Case 1: Chart spans International Date Line or Greenwich, Longitude
2053 // min/max is non-obvious.
2054 if ((m_LonMax * m_LonMin) < 0) // min/max are opposite signs
2055 {
2056 // Georeferencing is not yet available, so find the reference points
2057 // closest to min/max ply points
2058
2059 if (0 == nRefpoint) return false; // have to bail here
2060
2061 // for m_LonMax
2062 double min_dist_x = 360;
2063 int imaxclose = 0;
2064 for (int ic = 0; ic < nRefpoint; ic++) {
2065 double dist = sqrt(
2066 ((m_LatMax - pRefTable[ic].latr) * (m_LatMax - pRefTable[ic].latr)) +
2067 ((m_LonMax - pRefTable[ic].lonr) * (m_LonMax - pRefTable[ic].lonr)));
2068
2069 if (dist < min_dist_x) {
2070 min_dist_x = dist;
2071 imaxclose = ic;
2072 }
2073 }
2074
2075 // for m_LonMin
2076 double min_dist_n = 360;
2077 int iminclose = 0;
2078 for (int id = 0; id < nRefpoint; id++) {
2079 double dist = sqrt(
2080 ((m_LatMin - pRefTable[id].latr) * (m_LatMin - pRefTable[id].latr)) +
2081 ((m_LonMin - pRefTable[id].lonr) * (m_LonMin - pRefTable[id].lonr)));
2082
2083 if (dist < min_dist_n) {
2084 min_dist_n = dist;
2085 iminclose = id;
2086 }
2087 }
2088
2089 // Is this chart crossing IDL or Greenwich?
2090 // Make the check
2091 if (pRefTable[imaxclose].xr < pRefTable[iminclose].xr) {
2092 // This chart crosses IDL and needs a flip, meaning that all negative
2093 // longitudes need to be normalized and the min/max relcalculated This
2094 // code added to correct non-rectangular charts crossing IDL, such as
2095 // nz14605.kap
2096
2097 m_LonMax = -360.0;
2098 m_LonMin = 360.0;
2099 m_LatMax = -90.0;
2100 m_LatMin = 90.0;
2101
2102 Plypoint *ppp =
2103 (Plypoint *)GetCOVRTableHead(0); // Normalize the plypoints
2104 int cnPlypoint = GetCOVRTablenPoints(0);
2105
2106 for (int u = 0; u < cnPlypoint; u++) {
2107 if (ppp->lnp < 0.) ppp->lnp += 360.;
2108
2109 if (ppp->lnp > m_LonMax) m_LonMax = ppp->lnp;
2110 if (ppp->lnp < m_LonMin) m_LonMin = ppp->lnp;
2111
2112 if (ppp->ltp > m_LatMax) m_LatMax = ppp->ltp;
2113 if (ppp->ltp < m_LatMin) m_LatMin = ppp->ltp;
2114
2115 ppp++;
2116 }
2117 }
2118 }
2119
2120 // Case 2 Lons are both < -180, which means the extent will be reported
2121 // incorrectly and the plypoint structure will be wrong This case is seen
2122 // first on 81004_1.KAP, (Mariannas)
2123
2124 if ((m_LonMax < -180.) && (m_LonMin < -180.)) {
2125 m_LonMin += 360.; // Normalize the extents
2126 m_LonMax += 360.;
2127
2128 Plypoint *ppp = (Plypoint *)GetCOVRTableHead(0); // Normalize the plypoints
2129 int cnPlypoint = GetCOVRTablenPoints(0);
2130
2131 for (int u = 0; u < cnPlypoint; u++) {
2132 ppp->lnp += 360.;
2133 ppp++;
2134 }
2135 }
2136
2137 return true;
2138}
2139
2140void ChartBaseBSB::SetColorScheme(ColorScheme cs, bool bApplyImmediate) {
2141 // Here we convert (subjectively) the Global ColorScheme
2142 // to an appropriate BSB_Color_Capability index.
2143
2144 switch (cs) {
2145 case GLOBAL_COLOR_SCHEME_RGB:
2146 m_mapped_color_index = COLOR_RGB_DEFAULT;
2147 break;
2148 case GLOBAL_COLOR_SCHEME_DAY:
2149 m_mapped_color_index = DAY;
2150 break;
2151 case GLOBAL_COLOR_SCHEME_DUSK:
2152 m_mapped_color_index = DUSK;
2153 break;
2154 case GLOBAL_COLOR_SCHEME_NIGHT:
2155 m_mapped_color_index = NIGHT;
2156 break;
2157 default:
2158 m_mapped_color_index = DAY;
2159 break;
2160 }
2161
2162 pPalette = GetPalettePtr(m_mapped_color_index);
2163
2164 m_global_color_scheme = cs;
2165
2166 // Force a cache dump in a simple sideways manner
2167 if (bApplyImmediate) {
2168 m_cached_scale_ppm = 1.0;
2169 }
2170
2171 // Force a new thumbnail
2172 if (pThumbData) pThumbData->pDIBThumb = NULL;
2173}
2174
2175wxBitmap *ChartBaseBSB::CreateThumbnail(int tnx, int tny, ColorScheme cs) {
2176 // Calculate the size and divisors
2177
2178 int divx = wxMax(1, Size_X / (4 * tnx));
2179 int divy = wxMax(1, Size_Y / (4 * tny));
2180
2181 int div_factor = std::min(divx, divy);
2182
2183 int des_width = Size_X / div_factor;
2184 int des_height = Size_Y / div_factor;
2185
2186 wxRect gts;
2187 gts.x = 0; // full chart
2188 gts.y = 0;
2189 gts.width = Size_X;
2190 gts.height = Size_Y;
2191
2192 int this_bpp = 24; // for wxImage
2193 // Allocate the pixel storage needed for one line of chart bits
2194 unsigned char *pLineT = (unsigned char *)malloc((Size_X + 1) * BPP / 8);
2195
2196 // Scale the data quickly
2197 unsigned char *pPixTN =
2198 (unsigned char *)malloc(des_width * des_height * this_bpp / 8);
2199
2200 int ix = 0;
2201 int iy = 0;
2202 int iyd = 0;
2203 int ixd = 0;
2204 int yoffd;
2205 unsigned char *pxs;
2206 unsigned char *pxd;
2207
2208 // Temporarily set the color scheme
2209 ColorScheme cs_tmp = m_global_color_scheme;
2210 SetColorScheme(cs, false);
2211
2212 while (iyd < des_height) {
2213 if (0 == BSBGetScanline(pLineT, iy, 0, Size_X, 1)) // get a line
2214 {
2215 free(pLineT);
2216 free(pPixTN);
2217 return NULL;
2218 }
2219
2220 yoffd = iyd * des_width * this_bpp / 8; // destination y
2221
2222 ix = 0;
2223 ixd = 0;
2224 while (ixd < des_width) {
2225 pxs = pLineT + (ix * BPP / 8);
2226 pxd = pPixTN + (yoffd + (ixd * this_bpp / 8));
2227 *pxd++ = *pxs++;
2228 *pxd++ = *pxs++;
2229 *pxd = *pxs;
2230
2231 ix += div_factor;
2232 ixd++;
2233 }
2234
2235 iy += div_factor;
2236 iyd++;
2237 }
2238
2239 free(pLineT);
2240
2241 // Reset ColorScheme
2242 SetColorScheme(cs_tmp, false);
2243
2244 wxBitmap *retBMP;
2245
2246#ifdef ocpnUSE_ocpnBitmap
2247 wxBitmap *bmx2 = new ocpnBitmap(pPixTN, des_width, des_height, -1);
2248 wxImage imgx2 = bmx2->ConvertToImage();
2249 imgx2.Rescale(des_width / 4, des_height / 4, wxIMAGE_QUALITY_HIGH);
2250 retBMP = new wxBitmap(imgx2);
2251 delete bmx2;
2252#else
2253 wxImage thumb_image(des_width, des_height, pPixTN, true);
2254 thumb_image.Rescale(des_width / 4, des_height / 4, wxIMAGE_QUALITY_HIGH);
2255 retBMP = new wxBitmap(thumb_image);
2256#endif
2257
2258 free(pPixTN);
2259
2260 return retBMP;
2261}
2262
2263//-------------------------------------------------------------------------------------------------
2264// Get the Chart thumbnail data structure
2265// Creating the thumbnail bitmap as required
2266//-------------------------------------------------------------------------------------------------
2267
2268ThumbData *ChartBaseBSB::GetThumbData(int tnx, int tny, float lat, float lon) {
2269 // Create the bitmap if needed
2270 if (!pThumbData->pDIBThumb)
2271 pThumbData->pDIBThumb = CreateThumbnail(tnx, tny, m_global_color_scheme);
2272
2273 pThumbData->Thumb_Size_X = tnx;
2274 pThumbData->Thumb_Size_Y = tny;
2275
2276 // Plot the supplied Lat/Lon on the thumbnail
2277 int divx = Size_X / tnx;
2278 int divy = Size_Y / tny;
2279
2280 int div_factor = std::min(divx, divy);
2281
2282 double pixx, pixy;
2283
2284 // Using a temporary synthetic ViewPort and source rectangle,
2285 // calculate the ships position on the thumbnail
2286 ViewPort tvp;
2287 tvp.pix_width = tnx;
2288 tvp.pix_height = tny;
2289 tvp.view_scale_ppm = GetPPM() / div_factor;
2290 wxRect trex = Rsrc;
2291 Rsrc.x = 0;
2292 Rsrc.y = 0;
2293 latlong_to_pix_vp(lat, lon, pixx, pixy, tvp);
2294 Rsrc = trex;
2295
2296 pThumbData->ShipX = pixx; // / div_factor;
2297 pThumbData->ShipY = pixy; // / div_factor;
2298
2299 return pThumbData;
2300}
2301
2302bool ChartBaseBSB::UpdateThumbData(double lat, double lon) {
2303 // Plot the supplied Lat/Lon on the thumbnail
2304 // Return TRUE if the pixel location of ownship has changed
2305
2306 int divx = Size_X / pThumbData->Thumb_Size_X;
2307 int divy = Size_Y / pThumbData->Thumb_Size_Y;
2308
2309 int div_factor = std::min(divx, divy);
2310
2311 double pixx_test, pixy_test;
2312
2313 // Using a temporary synthetic ViewPort and source rectangle,
2314 // calculate the ships position on the thumbnail
2315 ViewPort tvp;
2316 tvp.pix_width = pThumbData->Thumb_Size_X;
2317 tvp.pix_height = pThumbData->Thumb_Size_Y;
2318 tvp.view_scale_ppm = GetPPM() / div_factor;
2319 wxRect trex = Rsrc;
2320 Rsrc.x = 0;
2321 Rsrc.y = 0;
2322 latlong_to_pix_vp(lat, lon, pixx_test, pixy_test, tvp);
2323 Rsrc = trex;
2324
2325 if ((pixx_test != pThumbData->ShipX) || (pixy_test != pThumbData->ShipY)) {
2326 pThumbData->ShipX = pixx_test;
2327 pThumbData->ShipY = pixy_test;
2328 return TRUE;
2329 } else
2330 return FALSE;
2331}
2332
2333//-----------------------------------------------------------------------
2334// Pixel to Lat/Long Conversion helpers
2335//-----------------------------------------------------------------------
2336static double polytrans(double *coeff, double lon, double lat);
2337
2338int ChartBaseBSB::vp_pix_to_latlong(ViewPort &vp, double pixx, double pixy,
2339 double *plat, double *plon) {
2340 if (bHaveEmbeddedGeoref) {
2341 double raster_scale = GetPPM() / vp.view_scale_ppm;
2342
2343 double px = pixx * raster_scale + Rsrc.x;
2344 double py = pixy * raster_scale + Rsrc.y;
2345 // pix_to_latlong(px, py, plat, plon);
2346
2347 if (1) {
2348 double lon = polytrans(pwx, px, py);
2349 lon = (lon < 0) ? lon + m_cph : lon - m_cph;
2350 *plon = lon - m_lon_datum_adjust;
2351 *plat = polytrans(pwy, px, py) - m_lat_datum_adjust;
2352 }
2353
2354 return 0;
2355 } else {
2356 double slat, slon;
2357 double xp, yp;
2358
2359 if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
2360 // Use Projected Polynomial algorithm
2361
2362 double raster_scale = GetPPM() / vp.view_scale_ppm;
2363
2364 // Apply poly solution to vp center point
2365 double easting, northing;
2366 toTM(vp.clat + m_lat_datum_adjust, vp.clon + m_lon_datum_adjust,
2367 m_proj_lat, m_proj_lon, &easting, &northing);
2368 double xc = polytrans(cPoints.wpx, easting, northing);
2369 double yc = polytrans(cPoints.wpy, easting, northing);
2370
2371 // convert screen pixels to chart pixmap relative
2372 double px = xc + (pixx - (vp.pix_width / 2)) * raster_scale;
2373 double py = yc + (pixy - (vp.pix_height / 2)) * raster_scale;
2374
2375 // Apply polynomial solution to chart relative pixels to get e/n
2376 double east = polytrans(cPoints.pwx, px, py);
2377 double north = polytrans(cPoints.pwy, px, py);
2378
2379 // Apply inverse Projection to get lat/lon
2380 double lat, lon;
2381 fromTM(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2382
2383 // Datum adjustments.....
2384 //?? lon = (lon < 0) ? lon + m_cph : lon - m_cph;
2385 double slon_p = lon - m_lon_datum_adjust;
2386 double slat_p = lat - m_lat_datum_adjust;
2387
2388 // printf("%8g %8g %8g %8g %g\n", slat, slat_p, slon,
2389 // slon_p, slon - slon_p);
2390 slon = slon_p;
2391 slat = slat_p;
2392
2393 } else if (m_projection == PROJECTION_MERCATOR) {
2394 // Use Projected Polynomial algorithm
2395
2396 double raster_scale = GetPPM() / vp.view_scale_ppm;
2397
2398 // Apply poly solution to vp center point
2399 double easting, northing;
2400 toSM_ECC(vp.clat + m_lat_datum_adjust, vp.clon + m_lon_datum_adjust,
2401 m_proj_lat, m_proj_lon, &easting, &northing);
2402 double xc = polytrans(cPoints.wpx, easting, northing);
2403 double yc = polytrans(cPoints.wpy, easting, northing);
2404
2405 // convert screen pixels to chart pixmap relative
2406 double px = xc + (pixx - (vp.pix_width / 2)) * raster_scale;
2407 double py = yc + (pixy - (vp.pix_height / 2)) * raster_scale;
2408
2409 // Apply polynomial solution to chart relative pixels to get e/n
2410 double east = polytrans(cPoints.pwx, px, py);
2411 double north = polytrans(cPoints.pwy, px, py);
2412
2413 // Apply inverse Projection to get lat/lon
2414 double lat, lon;
2415 fromSM_ECC(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2416
2417 // Make Datum adjustments.....
2418 double slon_p = lon - m_lon_datum_adjust;
2419 double slat_p = lat - m_lat_datum_adjust;
2420
2421 slon = slon_p;
2422 slat = slat_p;
2423
2424 // printf("vp.clon %g xc %g px %g east %g
2425 // \n", vp.clon, xc, px, east);
2426
2427 } else if (m_projection == PROJECTION_POLYCONIC) {
2428 // Use Projected Polynomial algorithm
2429
2430 double raster_scale = GetPPM() / vp.view_scale_ppm;
2431
2432 // Apply poly solution to vp center point
2433 double easting, northing;
2434 toPOLY(vp.clat + m_lat_datum_adjust, vp.clon + m_lon_datum_adjust,
2435 m_proj_lat, m_proj_lon, &easting, &northing);
2436 double xc = polytrans(cPoints.wpx, easting, northing);
2437 double yc = polytrans(cPoints.wpy, easting, northing);
2438
2439 // convert screen pixels to chart pixmap relative
2440 double px = xc + (pixx - (vp.pix_width / 2)) * raster_scale;
2441 double py = yc + (pixy - (vp.pix_height / 2)) * raster_scale;
2442
2443 // Apply polynomial solution to chart relative pixels to get e/n
2444 double east = polytrans(cPoints.pwx, px, py);
2445 double north = polytrans(cPoints.pwy, px, py);
2446
2447 // Apply inverse Projection to get lat/lon
2448 double lat, lon;
2449 fromPOLY(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2450
2451 // Make Datum adjustments.....
2452 double slon_p = lon - m_lon_datum_adjust;
2453 double slat_p = lat - m_lat_datum_adjust;
2454
2455 slon = slon_p;
2456 slat = slat_p;
2457
2458 } else {
2459 // Use a Mercator estimator, with Eccentricity corrrection applied
2460 int dx = pixx - (vp.pix_width / 2);
2461 int dy = (vp.pix_height / 2) - pixy;
2462
2463 xp = (dx * cos(vp.skew)) - (dy * sin(vp.skew));
2464 yp = (dy * cos(vp.skew)) + (dx * sin(vp.skew));
2465
2466 double d_east = xp / vp.view_scale_ppm;
2467 double d_north = yp / vp.view_scale_ppm;
2468
2469 fromSM_ECC(d_east, d_north, vp.clat, vp.clon, &slat, &slon);
2470 }
2471
2472 *plat = slat;
2473
2474 if (slon < -180.)
2475 slon += 360.;
2476 else if (slon > 180.)
2477 slon -= 360.;
2478 *plon = slon;
2479
2480 return 0;
2481 }
2482}
2483
2484int ChartBaseBSB::latlong_to_pix_vp(double lat, double lon, double &pixx,
2485 double &pixy, ViewPort &vp) {
2486 double alat, alon;
2487
2488 if (bHaveEmbeddedGeoref) {
2489 double alat, alon;
2490
2491 alon = lon + m_lon_datum_adjust;
2492 alat = lat + m_lat_datum_adjust;
2493
2494 AdjustLongitude(alon);
2495
2496 if (1) {
2497 /* change longitude phase (CPH) */
2498 double lonp = (alon < 0) ? alon + m_cph : alon - m_cph;
2499 double xd = polytrans(wpx, lonp, alat);
2500 double yd = polytrans(wpy, lonp, alat);
2501
2502 double raster_scale = GetPPM() / vp.view_scale_ppm;
2503
2504 pixx = (xd - Rsrc.x) / raster_scale;
2505 pixy = (yd - Rsrc.y) / raster_scale;
2506
2507 return 0;
2508 }
2509 } else {
2510 double easting, northing;
2511 double xlon = lon;
2512
2513 // Make sure lon and lon0 are same phase
2514 /*
2515 if((xlon * vp.clon) < 0.)
2516 {
2517 if(xlon < 0.)
2518 xlon += 360.;
2519 else
2520 xlon -= 360.;
2521 }
2522
2523 if(fabs(xlon - vp.clon) > 180.)
2524 {
2525 if(xlon > vp.clon)
2526 xlon -= 360.;
2527 else
2528 xlon += 360.;
2529 }
2530 */
2531
2532 if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
2533 // Use Projected Polynomial algorithm
2534
2535 alon = lon + m_lon_datum_adjust;
2536 alat = lat + m_lat_datum_adjust;
2537
2538 // Get e/n from TM Projection
2539 toTM(alat, alon, m_proj_lat, m_proj_lon, &easting, &northing);
2540
2541 // Apply poly solution to target point
2542 double xd = polytrans(cPoints.wpx, easting, northing);
2543 double yd = polytrans(cPoints.wpy, easting, northing);
2544
2545 // Apply poly solution to vp center point
2546 toTM(vp.clat + m_lat_datum_adjust, vp.clon + m_lon_datum_adjust,
2547 m_proj_lat, m_proj_lon, &easting, &northing);
2548 double xc = polytrans(cPoints.wpx, easting, northing);
2549 double yc = polytrans(cPoints.wpy, easting, northing);
2550
2551 // Calculate target point relative to vp center
2552 double raster_scale = GetPPM() / vp.view_scale_ppm;
2553
2554 double xs = xc - vp.pix_width * raster_scale / 2;
2555 double ys = yc - vp.pix_height * raster_scale / 2;
2556
2557 pixx = (xd - xs) / raster_scale;
2558 pixy = (yd - ys) / raster_scale;
2559
2560 } else if (m_projection == PROJECTION_MERCATOR) {
2561 // Use Projected Polynomial algorithm
2562
2563 alon = lon + m_lon_datum_adjust;
2564 alat = lat + m_lat_datum_adjust;
2565
2566 // Get e/n from Projection
2567 xlon = alon;
2568 AdjustLongitude(xlon);
2569 toSM_ECC(alat, xlon, m_proj_lat, m_proj_lon, &easting, &northing);
2570
2571 // Apply poly solution to target point
2572 double xd = polytrans(cPoints.wpx, easting, northing);
2573 double yd = polytrans(cPoints.wpy, easting, northing);
2574
2575 // Apply poly solution to vp center point
2576 double xlonc = vp.clon;
2577 AdjustLongitude(xlonc);
2578
2579 toSM_ECC(vp.clat + m_lat_datum_adjust, xlonc + m_lon_datum_adjust,
2580 m_proj_lat, m_proj_lon, &easting, &northing);
2581 double xc = polytrans(cPoints.wpx, easting, northing);
2582 double yc = polytrans(cPoints.wpy, easting, northing);
2583
2584 // Calculate target point relative to vp center
2585 double raster_scale = GetPPM() / vp.view_scale_ppm;
2586
2587 double xs = xc - vp.pix_width * raster_scale / 2;
2588 double ys = yc - vp.pix_height * raster_scale / 2;
2589
2590 pixx = (xd - xs) / raster_scale;
2591 pixy = (yd - ys) / raster_scale;
2592
2593 } else if (m_projection == PROJECTION_POLYCONIC) {
2594 // Use Projected Polynomial algorithm
2595
2596 alon = lon + m_lon_datum_adjust;
2597 alat = lat + m_lat_datum_adjust;
2598
2599 // Get e/n from Projection
2600 xlon = AdjustLongitude(alon);
2601 toPOLY(alat, xlon, m_proj_lat, m_proj_lon, &easting, &northing);
2602
2603 // Apply poly solution to target point
2604 double xd = polytrans(cPoints.wpx, easting, northing);
2605 double yd = polytrans(cPoints.wpy, easting, northing);
2606
2607 // Apply poly solution to vp center point
2608 double xlonc = AdjustLongitude(vp.clon);
2609
2610 toPOLY(vp.clat + m_lat_datum_adjust, xlonc + m_lon_datum_adjust,
2611 m_proj_lat, m_proj_lon, &easting, &northing);
2612 double xc = polytrans(cPoints.wpx, easting, northing);
2613 double yc = polytrans(cPoints.wpy, easting, northing);
2614
2615 // Calculate target point relative to vp center
2616 double raster_scale = GetPPM() / vp.view_scale_ppm;
2617
2618 double xs = xc - vp.pix_width * raster_scale / 2;
2619 double ys = yc - vp.pix_height * raster_scale / 2;
2620
2621 pixx = (xd - xs) / raster_scale;
2622 pixy = (yd - ys) / raster_scale;
2623
2624 } else {
2625 toSM_ECC(lat, xlon, vp.clat, vp.clon, &easting, &northing);
2626
2627 double epix = easting * vp.view_scale_ppm;
2628 double npix = northing * vp.view_scale_ppm;
2629
2630 double dx = epix * cos(vp.skew) + npix * sin(vp.skew);
2631 double dy = npix * cos(vp.skew) - epix * sin(vp.skew);
2632
2633 pixx = ((double)vp.pix_width / 2) + dx;
2634 pixy = ((double)vp.pix_height / 2) - dy;
2635 }
2636 return 0;
2637 }
2638
2639 return 1;
2640}
2641
2642void ChartBaseBSB::latlong_to_chartpix(double lat, double lon, double &pixx,
2643 double &pixy) {
2644 double alat, alon;
2645 pixx = 0.0;
2646 pixy = 0.0;
2647
2648 if (bHaveEmbeddedGeoref) {
2649 double alat, alon;
2650
2651 alon = lon + m_lon_datum_adjust;
2652 alat = lat + m_lat_datum_adjust;
2653
2654 alon = AdjustLongitude(alon);
2655
2656 /* change longitude phase (CPH) */
2657 double lonp = (alon < 0) ? alon + m_cph : alon - m_cph;
2658 pixx = polytrans(wpx, lonp, alat);
2659 pixy = polytrans(wpy, lonp, alat);
2660 } else {
2661 double easting, northing;
2662 double xlon = lon;
2663
2664 if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
2665 // Use Projected Polynomial algorithm
2666
2667 alon = lon + m_lon_datum_adjust;
2668 alat = lat + m_lat_datum_adjust;
2669
2670 // Get e/n from TM Projection
2671 toTM(alat, alon, m_proj_lat, m_proj_lon, &easting, &northing);
2672
2673 // Apply poly solution to target point
2674 pixx = polytrans(cPoints.wpx, easting, northing);
2675 pixy = polytrans(cPoints.wpy, easting, northing);
2676
2677 } else if (m_projection == PROJECTION_MERCATOR) {
2678 // Use Projected Polynomial algorithm
2679
2680 alon = lon + m_lon_datum_adjust;
2681 alat = lat + m_lat_datum_adjust;
2682
2683 // Get e/n from Projection
2684 xlon = AdjustLongitude(alon);
2685
2686 toSM_ECC(alat, xlon, m_proj_lat, m_proj_lon, &easting, &northing);
2687
2688 // Apply poly solution to target point
2689 pixx = polytrans(cPoints.wpx, easting, northing);
2690 pixy = polytrans(cPoints.wpy, easting, northing);
2691
2692 } else if (m_projection == PROJECTION_POLYCONIC) {
2693 // Use Projected Polynomial algorithm
2694
2695 alon = lon + m_lon_datum_adjust;
2696 alat = lat + m_lat_datum_adjust;
2697
2698 // Get e/n from Projection
2699 xlon = AdjustLongitude(alon);
2700 toPOLY(alat, xlon, m_proj_lat, m_proj_lon, &easting, &northing);
2701
2702 // Apply poly solution to target point
2703 pixx = polytrans(cPoints.wpx, easting, northing);
2704 pixy = polytrans(cPoints.wpy, easting, northing);
2705 }
2706 }
2707}
2708
2709void ChartBaseBSB::chartpix_to_latlong(double pixx, double pixy, double *plat,
2710 double *plon) {
2711 if (bHaveEmbeddedGeoref) {
2712 double lon = polytrans(pwx, pixx, pixy);
2713 lon = (lon < 0) ? lon + m_cph : lon - m_cph;
2714 *plon = lon - m_lon_datum_adjust;
2715 *plat = polytrans(pwy, pixx, pixy) - m_lat_datum_adjust;
2716 } else {
2717 double slat, slon;
2718 if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
2719 // Use Projected Polynomial algorithm
2720
2721 // Apply polynomial solution to chart relative pixels to get e/n
2722 double east = polytrans(cPoints.pwx, pixx, pixy);
2723 double north = polytrans(cPoints.pwy, pixx, pixy);
2724
2725 // Apply inverse Projection to get lat/lon
2726 double lat, lon;
2727 fromTM(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2728
2729 // Datum adjustments.....
2730 //?? lon = (lon < 0) ? lon + m_cph : lon - m_cph;
2731 slon = lon - m_lon_datum_adjust;
2732 slat = lat - m_lat_datum_adjust;
2733
2734 } else if (m_projection == PROJECTION_MERCATOR) {
2735 // Use Projected Polynomial algorithm
2736 // Apply polynomial solution to chart relative pixels to get e/n
2737 double east = polytrans(cPoints.pwx, pixx, pixy);
2738 double north = polytrans(cPoints.pwy, pixx, pixy);
2739
2740 // Apply inverse Projection to get lat/lon
2741 double lat, lon;
2742 fromSM_ECC(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2743
2744 // Make Datum adjustments.....
2745 slon = lon - m_lon_datum_adjust;
2746 slat = lat - m_lat_datum_adjust;
2747 } else if (m_projection == PROJECTION_POLYCONIC) {
2748 // Use Projected Polynomial algorithm
2749 // Apply polynomial solution to chart relative pixels to get e/n
2750 double east = polytrans(cPoints.pwx, pixx, pixy);
2751 double north = polytrans(cPoints.pwy, pixx, pixy);
2752
2753 // Apply inverse Projection to get lat/lon
2754 double lat, lon;
2755 fromPOLY(east, north, m_proj_lat, m_proj_lon, &lat, &lon);
2756
2757 // Make Datum adjustments.....
2758 slon = lon - m_lon_datum_adjust;
2759 slat = lat - m_lat_datum_adjust;
2760
2761 } else {
2762 slon = 0.;
2763 slat = 0.;
2764 }
2765
2766 *plat = slat;
2767
2768 if (slon < -180.)
2769 slon += 360.;
2770 else if (slon > 180.)
2771 slon -= 360.;
2772 *plon = slon;
2773 }
2774}
2775
2776void ChartBaseBSB::ComputeSourceRectangle(const ViewPort &vp,
2777 wxRect *pSourceRect) {
2778 m_raster_scale_factor = GetRasterScaleFactor(vp);
2779 double xd, yd;
2780 latlong_to_chartpix(vp.clat, vp.clon, xd, yd);
2781
2782 wxRealPoint pos, size;
2783
2784 pos.x = xd - (vp.pix_width * m_raster_scale_factor / 2);
2785 pos.y = yd - (vp.pix_height * m_raster_scale_factor / 2);
2786
2787 size.x = vp.pix_width * m_raster_scale_factor;
2788 size.y = vp.pix_height * m_raster_scale_factor;
2789
2790 *pSourceRect =
2791 wxRect(wxRound(pos.x), wxRound(pos.y), wxRound(size.x), wxRound(size.y));
2792}
2793
2794double ChartBaseBSB::GetRasterScaleFactor(const ViewPort &vp) {
2795 // This funny contortion is necessary to allow scale factors < 1, i.e.
2796 // overzoom
2797 return (wxRound(100000 * GetPPM() / vp.view_scale_ppm)) / 100000.;
2798}
2799
2800void ChartBaseBSB::SetVPRasterParms(const ViewPort &vpt) {
2801 // Calculate the potential datum offset parameters for this viewport, if
2802 // not WGS84
2803
2804 if (m_datum_index == DATUM_INDEX_WGS84 ||
2805 m_datum_index == DATUM_INDEX_UNKNOWN) {
2806 m_lon_datum_adjust = (-m_dtm_lon) / 3600.;
2807 m_lat_datum_adjust = (-m_dtm_lat) / 3600.;
2808 }
2809
2810 else {
2811 double to_lat, to_lon;
2812 MolodenskyTransform(vpt.clat, vpt.clon, &to_lat, &to_lon, m_datum_index,
2813 DATUM_INDEX_WGS84);
2814 m_lon_datum_adjust = -(to_lon - vpt.clon);
2815 m_lat_datum_adjust = -(to_lat - vpt.clat);
2816 if (m_b_apply_dtm) {
2817 m_lon_datum_adjust -= m_dtm_lon / 3600.;
2818 m_lat_datum_adjust -= m_dtm_lat / 3600.;
2819 }
2820 }
2821
2822 ComputeSourceRectangle(vpt, &Rsrc);
2823
2824 if (vpt.IsValid()) m_vp_render_last = vpt;
2825}
2826
2827bool ChartBaseBSB::AdjustVP(ViewPort &vp_last, ViewPort &vp_proposed) {
2828 bool bInside = G_FloatPtInPolygon((MyFlPoint *)GetCOVRTableHead(0),
2829 GetCOVRTablenPoints(0), vp_proposed.clon,
2830 vp_proposed.clat);
2831 if (!bInside) return false;
2832
2833 int ret_val = 0;
2834 double binary_scale_factor = GetPPM() / vp_proposed.view_scale_ppm;
2835
2836 if (vp_last.IsValid()) {
2837 // We only need to adjust the VP if the cache is valid and potentially
2838 // usable, i.e. the scale factor is integer... The objective here is to
2839 // ensure that the VP center falls on an exact pixel boundary within the
2840 // cache
2841
2842 if (cached_image_ok && (binary_scale_factor > 1.0) &&
2843 (fabs(binary_scale_factor - wxRound(binary_scale_factor)) < 1e-5)) {
2844 if (m_b_cdebug)
2845 printf(" Possible Adjust VP for integer scale: %g\n",
2846 binary_scale_factor);
2847
2848 wxRect rprop;
2849 ComputeSourceRectangle(vp_proposed, &rprop);
2850
2851 double pixx, pixy;
2852 double lon_adj, lat_adj;
2853 latlong_to_pix_vp(vp_proposed.clat, vp_proposed.clon, pixx, pixy,
2854 vp_proposed);
2855 vp_pix_to_latlong(vp_proposed, pixx, pixy, &lat_adj, &lon_adj);
2856
2857 vp_proposed.clat = lat_adj;
2858 vp_proposed.clon = lon_adj;
2859 ret_val = 1;
2860 }
2861 }
2862
2863 return (ret_val > 0);
2864}
2865
2866bool ChartBaseBSB::IsRenderCacheable(wxRect &source, wxRect &dest) {
2867 double scale_x = (double)source.width / (double)dest.width;
2868
2869 if (scale_x <= 1.0) // overzoom
2870 {
2871 // if(m_b_cdebug)printf(" MISS<<<>>>GVUC: Overzoom\n");
2872 return false;
2873 }
2874
2875 // Using the cache only works for pure binary scale factors......
2876 if ((fabs(scale_x - wxRound(scale_x))) > .0001) {
2877 // if(m_b_cdebug)printf(" MISS<<<>>>GVUC: Not digital scale
2878 // test 1\n");
2879 return false;
2880 }
2881
2882 // Scale must be exactly digital...
2883 if ((int)(source.width / dest.width) != (int)wxRound(scale_x)) {
2884 // if(m_b_cdebug)printf(" MISS<<<>>>GVUC: Not digital scale
2885 // test 2\n");
2886 return false;
2887 }
2888
2889 return true;
2890}
2891
2892void ChartBaseBSB::GetValidCanvasRegion(const ViewPort &VPoint,
2893 OCPNRegion *pValidRegion) {
2894 SetVPRasterParms(VPoint);
2895
2896 double raster_scale = VPoint.view_scale_ppm / GetPPM();
2897
2898 int rxl, rxr;
2899 int ryb, ryt;
2900
2901 rxl = wxMax(-Rsrc.x * raster_scale, VPoint.rv_rect.x);
2902 rxr = wxMin((Size_X - Rsrc.x) * raster_scale,
2903 VPoint.rv_rect.width + VPoint.rv_rect.x);
2904
2905 ryt = wxMax(-Rsrc.y * raster_scale, VPoint.rv_rect.y);
2906 ryb = wxMin((Size_Y - Rsrc.y) * raster_scale,
2907 VPoint.rv_rect.height + VPoint.rv_rect.y);
2908
2909 pValidRegion->Clear();
2910 pValidRegion->Union(rxl, ryt, rxr - rxl, ryb - ryt);
2911}
2912
2913LLRegion ChartBaseBSB::GetValidRegion() {
2914 // should we cache this?
2915 double ll[8];
2916 chartpix_to_latlong(0, 0, ll + 0, ll + 1);
2917 chartpix_to_latlong(0, Size_Y, ll + 2, ll + 3);
2918 chartpix_to_latlong(Size_X, Size_Y, ll + 4, ll + 5);
2919 chartpix_to_latlong(Size_X, 0, ll + 6, ll + 7);
2920
2921 // for now don't allow raster charts to cross both 0 meridian and IDL
2922 // (complicated to deal with)
2923 for (int i = 1; i < 6; i += 2)
2924 if (fabs(ll[i] - ll[i + 2]) > 180) {
2925 // we detect crossing idl here, make all longitudes positive
2926 for (int i = 1; i < 8; i += 2)
2927 if (ll[i] < 0) ll[i] += 360;
2928 break;
2929 }
2930
2931 return LLRegion(4, ll);
2932}
2933
2934bool ChartBaseBSB::GetViewUsingCache(wxRect &source, wxRect &dest,
2935 const OCPNRegion &Region,
2936 ScaleTypeEnum scale_type) {
2937 wxRect s1;
2938 ScaleTypeEnum scale_type_corrected;
2939
2940 if (m_b_cdebug) printf(" source: %d %d\n", source.x, source.y);
2941 if (m_b_cdebug) printf(" cache: %d %d\n", cache_rect.x, cache_rect.y);
2942
2943 // Anything to do?
2944 if ((source == cache_rect) /*&& (cache_scale_method == scale_type)*/ &&
2945 (cached_image_ok)) {
2946 if (m_b_cdebug) printf(" GVUC: Cache is good, nothing to do\n");
2947 return false;
2948 }
2949
2950 double scale_x = (double)source.width / (double)dest.width;
2951
2952 if (m_b_cdebug) printf("GVUC: scale_x: %g\n", scale_x);
2953
2954 // Enforce a limit on bilinear scaling, for performance reasons
2955 scale_type_corrected = scale_type; // RENDER_LODEF; //scale_type;
2956 if (scale_x > m_bilinear_limit) scale_type_corrected = RENDER_LODEF;
2957
2958 {
2959 // if(b_cdebug)printf(" MISS<<<>>>GVUC: Intentional out\n");
2960 // return GetView( source, dest, scale_type_corrected );
2961 }
2962
2963 // Using the cache only works for pure binary scale factors......
2964 if ((fabs(scale_x - wxRound(scale_x))) > .0001) {
2965 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: Not digital scale test 1\n");
2966 return GetView(source, dest, scale_type_corrected);
2967 }
2968
2969 // scale_type_corrected = RENDER_LODEF;
2970
2971 if (!cached_image_ok) {
2972 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: Cache NOk\n");
2973 return GetView(source, dest, scale_type_corrected);
2974 }
2975
2976 if (scale_x < 1.0) // overzoom
2977 {
2978 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: Overzoom\n");
2979 return GetView(source, dest, scale_type_corrected);
2980 }
2981
2982 // Scale must be exactly digital...
2983 if ((int)(source.width / dest.width) != (int)wxRound(scale_x)) {
2984 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: Not digital scale test 2\n");
2985 return GetView(source, dest, scale_type_corrected);
2986 }
2987
2988 // Calculate the digital scale, e.g. 1,2,4,8,,,
2989 int cs1d = source.width / dest.width;
2990 if (abs(source.x - cache_rect.x) % cs1d) {
2991 if (m_b_cdebug)
2992 printf(" source.x: %d cache_rect.x: %d cs1d: %d\n", source.x,
2993 cache_rect.x, cs1d);
2994 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: x mismatch\n");
2995 return GetView(source, dest, scale_type_corrected);
2996 }
2997 if (abs(source.y - cache_rect.y) % cs1d) {
2998 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: y mismatch\n");
2999 return GetView(source, dest, scale_type_corrected);
3000 }
3001
3002 if (pPixCache && ((pPixCache->GetWidth() != dest.width) ||
3003 (pPixCache->GetHeight() != dest.height))) {
3004 if (m_b_cdebug) printf(" MISS<<<>>>GVUC: dest size mismatch\n");
3005 return GetView(source, dest, scale_type_corrected);
3006 }
3007
3008 int stride_rows =
3009 (source.y + source.height) - (cache_rect.y + cache_rect.height);
3010 int scaled_stride_rows = (int)(stride_rows / scale_x);
3011
3012 if (abs(stride_rows) >= source.height) // Pan more than one screen
3013 return GetView(source, dest, scale_type_corrected);
3014
3015 int stride_pixels =
3016 (source.x + source.width) - (cache_rect.x + cache_rect.width);
3017 int scaled_stride_pixels = (int)(stride_pixels / scale_x);
3018
3019 if (abs(stride_pixels) >= source.width) // Pan more than one screen
3020 return GetView(source, dest, scale_type_corrected);
3021
3022 if (m_b_cdebug) printf(" GVUC Using raster data cache\n");
3023
3024 ScaleTypeEnum pan_scale_type_x = scale_type_corrected;
3025 ScaleTypeEnum pan_scale_type_y = scale_type_corrected;
3026
3027 // "Blit" the valid pixels out of the way
3028 if (pPixCache) {
3029 int height = pPixCache->GetHeight();
3030 int width = pPixCache->GetWidth();
3031 int buffer_stride_bytes = pPixCache->GetLinePitch();
3032
3033 unsigned char *ps;
3034 unsigned char *pd;
3035
3036 if (stride_rows > 0) // pan down
3037 {
3038 ps = pPixCache->GetpData() +
3039 (abs(scaled_stride_rows) * buffer_stride_bytes);
3040 if (stride_pixels > 0) ps += scaled_stride_pixels * BPP / 8;
3041
3042 pd = pPixCache->GetpData();
3043 if (stride_pixels <= 0) pd += abs(scaled_stride_pixels) * BPP / 8;
3044
3045 for (int iy = 0; iy < (height - abs(scaled_stride_rows)); iy++) {
3046 memmove(pd, ps, (width - abs(scaled_stride_pixels)) * BPP / 8);
3047 ps += buffer_stride_bytes;
3048 pd += buffer_stride_bytes;
3049 }
3050
3051 } else {
3052 ps = pPixCache->GetpData() +
3053 ((height - abs(scaled_stride_rows) - 1) * buffer_stride_bytes);
3054 if (stride_pixels > 0) // make a hole on right
3055 ps += scaled_stride_pixels * BPP / 8;
3056
3057 pd = pPixCache->GetpData() + ((height - 1) * buffer_stride_bytes);
3058 if (stride_pixels <= 0) // make a hole on the left
3059 pd += abs(scaled_stride_pixels) * BPP / 8;
3060
3061 for (int iy = 0; iy < (height - abs(scaled_stride_rows)); iy++) {
3062 memmove(pd, ps, (width - abs(scaled_stride_pixels)) * BPP / 8);
3063 ps -= buffer_stride_bytes;
3064 pd -= buffer_stride_bytes;
3065 }
3066 }
3067
3068 // Y Pan
3069 if (source.y != cache_rect.y) {
3070 wxRect sub_dest = dest;
3071 sub_dest.height = abs(scaled_stride_rows);
3072
3073 if (stride_rows > 0) // pan down
3074 {
3075 sub_dest.y = height - scaled_stride_rows;
3076
3077 } else {
3078 sub_dest.y = 0;
3079 }
3080
3081 // Get the new bits needed
3082
3083 // A little optimization...
3084 // No sense in fetching bits that are not part of the ultimate render
3085 // region
3086 wxRegionContain rc = Region.Contains(sub_dest);
3087 if ((wxPartRegion == rc) || (wxInRegion == rc)) {
3088 GetAndScaleData(pPixCache->GetpData(), pPixCache->GetLength(), source,
3089 source.width, sub_dest, width, cs1d, pan_scale_type_y);
3090 }
3091 pPixCache->Update();
3092
3093 // Update the cached parameters, Y only
3094
3095 cache_rect.y = source.y;
3096 // cache_rect = source;
3097 cache_rect_scaled = dest;
3098 cached_image_ok = 1;
3099
3100 } // Y Pan
3101
3102 // X Pan
3103 if (source.x != cache_rect.x) {
3104 wxRect sub_dest = dest;
3105 sub_dest.width = abs(scaled_stride_pixels);
3106
3107 if (stride_pixels > 0) // pan right
3108 {
3109 sub_dest.x = width - scaled_stride_pixels;
3110 } else // pan left
3111 {
3112 sub_dest.x = 0;
3113 }
3114
3115 // Get the new bits needed
3116
3117 // A little optimization...
3118 // No sense in fetching bits that are not part of the ultimate render
3119 // region
3120 wxRegionContain rc = Region.Contains(sub_dest);
3121 if ((wxPartRegion == rc) || (wxInRegion == rc)) {
3122 GetAndScaleData(pPixCache->GetpData(), pPixCache->GetLength(), source,
3123 source.width, sub_dest, width, cs1d, pan_scale_type_x);
3124 }
3125
3126 pPixCache->Update();
3127
3128 // Update the cached parameters
3129 cache_rect = source;
3130 cache_rect_scaled = dest;
3131 cached_image_ok = 1;
3132
3133 } // X pan
3134
3135 return true;
3136 }
3137 return false;
3138}
3139
3140int s_dc;
3141
3142bool ChartBaseBSB::RenderViewOnDC(wxMemoryDC &dc, const ViewPort &VPoint) {
3143 SetVPRasterParms(VPoint);
3144
3145 OCPNRegion rgn(0, 0, VPoint.pix_width, VPoint.pix_height);
3146
3147 bool bsame_region = (rgn == m_last_region); // only want to do this once
3148
3149 if (!bsame_region) cached_image_ok = false;
3150
3151 m_last_region = rgn;
3152
3153 return RenderRegionViewOnDC(dc, VPoint, rgn);
3154}
3155
3156bool ChartBaseBSB::RenderRegionViewOnGL(const wxGLContext &glc,
3157 const ViewPort &VPoint,
3158 const OCPNRegion &RectRegion,
3159 const LLRegion &Region) {
3160 return true;
3161}
3162
3163bool ChartBaseBSB::RenderRegionViewOnDC(wxMemoryDC &dc, const ViewPort &VPoint,
3164 const OCPNRegion &Region) {
3165 SetVPRasterParms(VPoint);
3166
3167 wxRect dest(0, 0, VPoint.pix_width, VPoint.pix_height);
3168 // double factor = ((double)Rsrc.width)/((double)dest.width);
3169 double factor = GetRasterScaleFactor(VPoint);
3170 if (m_b_cdebug) {
3171 printf("%d RenderRegion ScaleType: %d factor: %g\n", s_dc++,
3172 RENDER_HIDEF, factor);
3173 printf("Rect list:\n");
3174 OCPNRegionIterator upd(Region); // get the requested rect list
3175 while (upd.HaveRects()) {
3176 wxRect rect = upd.GetRect();
3177 printf(" %d %d %d %d\n", rect.x, rect.y, rect.width, rect.height);
3178 upd.NextRect();
3179 }
3180 }
3181
3182 // Invalidate the cache if the scale has changed or the viewport size has
3183 // changed....
3184 if ((fabs(m_cached_scale_ppm - VPoint.view_scale_ppm) > 1e-9) ||
3185 (m_last_vprect != dest)) {
3186 cached_image_ok = false;
3187 m_vp_render_last.Invalidate();
3188 }
3189 /*
3190 if(pPixCache)
3191 {
3192 if((pPixCache->GetWidth() != dest.width) ||
3193 (pPixCache->GetHeight() != dest.height))
3194 {
3195 delete pPixCache;
3196 pPixCache = new PixelCache(dest.width, dest.height, BPP);
3197 }
3198 }
3199 else
3200 pPixCache = new PixelCache(dest.width, dest.height, BPP);
3201 */
3202
3203 m_cached_scale_ppm = VPoint.view_scale_ppm;
3204 m_last_vprect = dest;
3205
3206 if (cached_image_ok) {
3207 // Anything to do?
3208 bool bsame_region = (Region == m_last_region); // only want to do this once
3209
3210 if ((bsame_region) && (Rsrc == cache_rect)) {
3211 pPixCache->SelectIntoDC(dc);
3212 if (m_b_cdebug) printf(" Using Current PixelCache\n");
3213 return false;
3214 }
3215 }
3216
3217 m_last_region = Region;
3218
3219 // Analyze the region requested
3220 // When rendering complex regions, (more than say 4 rectangles)
3221 // .OR. small proportions, then rectangle rendering may be faster
3222 // Also true if the scale is less than near unity, or overzoom.
3223 // This will be the case for backgrounds of the quilt.
3224
3225 /* Update for Version 2.4.0
3226 This logic seems flawed, at least for quilts which contain charts having
3227 non-rectangular coverage areas. These quilt regions decompose to ...LOTS... of
3228 rectangles, most of which are 1 pixel in height. This is very slow, due to the
3229 overhead of GetAndScaleData(). However, remember that overzoom never uses the
3230 cache, nor does non-binary scale factors.. So, we check to see if this is a
3231 cacheable render, and that the number of rectangles is "reasonable"
3232 */
3233
3234 // Get the region rectangle count
3235
3236 int n_rect = 0;
3237 OCPNRegionIterator upd(Region); // get the requested rect list
3238 while (upd.HaveRects()) {
3239 n_rect++;
3240 upd.NextRect();
3241 }
3242
3243 if ((!IsRenderCacheable(Rsrc, dest) && (n_rect > 4) && (n_rect < 20)) ||
3244 (factor < 1)) {
3245 if (m_b_cdebug)
3246 printf(" RenderRegion by rect iterator n_rect: %d\n", n_rect);
3247
3248 // Verify that the persistent pixel cache is at least as large as the
3249 // largest rectangle in the region
3250 wxRect dest_check_rect = dest;
3251 OCPNRegionIterator upd_check(Region); // get the requested rect list
3252 while (upd_check.HaveRects()) {
3253 wxRect rect = upd_check.GetRect();
3254 dest_check_rect.Union(rect);
3255 upd_check.NextRect();
3256 }
3257
3258 if (pPixCache) {
3259 if ((pPixCache->GetWidth() != dest_check_rect.width) ||
3260 (pPixCache->GetHeight() != dest_check_rect.height)) {
3261 delete pPixCache;
3262 pPixCache =
3263 new PixelCache(dest_check_rect.width, dest_check_rect.height, BPP);
3264 }
3265 } else
3266 pPixCache =
3267 new PixelCache(dest_check_rect.width, dest_check_rect.height, BPP);
3268
3269 ScaleTypeEnum ren_type = RENDER_LODEF;
3270
3271 // Decompose the region into rectangles, and fetch them into the target
3272 // dc
3273 OCPNRegionIterator upd(Region); // get the requested rect list
3274 int ir = 0;
3275 while (upd.HaveRects()) {
3276 wxRect rect = upd.GetRect();
3277
3278 // Floating point math can lead to negative rectangle origin.
3279 // If this happens, we arbitrarily shift the rectangle to be positive
3280 // semidefinite. This will cause at most a 1 pixlel error onscreen.
3281 if (rect.y < 0) rect.Offset(0, -rect.y);
3282 if (rect.x < 0) rect.Offset(-rect.x, 0);
3283
3284 GetAndScaleData(pPixCache->GetpData(), pPixCache->GetLength(), Rsrc,
3285 Rsrc.width, rect, pPixCache->GetWidth(), factor,
3286 ren_type);
3287
3288 ir++;
3289 upd.NextRect();
3290 ;
3291 }
3292
3293 pPixCache->Update();
3294
3295 // Update cache parameters
3296 cache_rect = Rsrc;
3297 cache_scale_method = ren_type;
3298 cached_image_ok = false; // Never cache this type of render
3299
3300 // Select the data into the dc
3301 pPixCache->SelectIntoDC(dc);
3302
3303 return true;
3304 }
3305
3306 // Default is to try using the cache
3307
3308 if (pPixCache) {
3309 if ((pPixCache->GetWidth() != dest.width) ||
3310 (pPixCache->GetHeight() != dest.height)) {
3311 delete pPixCache;
3312 pPixCache = new PixelCache(dest.width, dest.height, BPP);
3313 }
3314 } else
3315 pPixCache = new PixelCache(dest.width, dest.height, BPP);
3316
3317 if (m_b_cdebug) printf(" Render Region By GVUC\n");
3318
3319 // A performance enhancement.....
3320 ScaleTypeEnum scale_type_zoom = RENDER_HIDEF;
3321 double binary_scale_factor = VPoint.view_scale_ppm / GetPPM();
3322 if (binary_scale_factor < .20) scale_type_zoom = RENDER_LODEF;
3323
3324 bool bnewview = GetViewUsingCache(Rsrc, dest, Region, scale_type_zoom);
3325
3326 // Select the data into the dc
3327 pPixCache->SelectIntoDC(dc);
3328
3329 return bnewview;
3330}
3331
3332wxImage *ChartBaseBSB::GetImage() {
3333 int img_size_x = ((Size_X >> 2) * 4) + 4;
3334 wxImage *img = new wxImage(img_size_x, Size_Y, false);
3335
3336 unsigned char *ppnx = img->GetData();
3337
3338 for (int i = 0; i < Size_Y; i++) {
3339 wxRect source_rect(0, i, Size_X, 1);
3340 wxRect dest_rect(0, 0, Size_X, 1);
3341
3342 GetAndScaleData(img->GetData(), img_size_x * Size_Y * 3, source_rect,
3343 Size_X, dest_rect, Size_X, 1.0, RENDER_HIDEF);
3344
3345 ppnx += img_size_x * 3;
3346 }
3347
3348 return img;
3349}
3350
3351bool ChartBaseBSB::GetView(wxRect &source, wxRect &dest,
3352 ScaleTypeEnum scale_type) {
3353 assert(pPixCache != 0);
3354 // PixelCache *pPixCacheTemp = new PixelCache(dest.width, dest.height,
3355 // BPP);
3356
3357 // Get and Rescale the data directly into the temporary PixelCache data
3358 // buffer
3359 double factor = ((double)source.width) / ((double)dest.width);
3360
3361 /*
3362 if(!GetAndScaleData(&ppnx, source, source.width, dest, dest.width,
3363 factor, scale_type))
3364 {
3365 delete pPixCacheTemp; // Some error, retain
3366 old cache return false;
3367 }
3368 else
3369 {
3370 delete pPixCache; // new cache is OK
3371 pPixCache = pPixCacheTemp;
3372 }
3373 */
3374 GetAndScaleData(pPixCache->GetpData(), pPixCache->GetLength(), source,
3375 source.width, dest, dest.width, factor, scale_type);
3376 pPixCache->Update();
3377
3378 // Update cache parameters
3379
3380 cache_rect = source;
3381 cache_rect_scaled = dest;
3382 cache_scale_method = scale_type;
3383
3384 cached_image_ok = 1;
3385
3386 return TRUE;
3387}
3388
3389bool ChartBaseBSB::GetAndScaleData(unsigned char *ppn, size_t data_size,
3390 wxRect &source, int source_stride,
3391 wxRect &dest, int dest_stride,
3392 double scale_factor,
3393 ScaleTypeEnum scale_type) {
3394 unsigned char *s_data = NULL;
3395
3396 double factor = scale_factor;
3397 int Factor = (int)factor;
3398
3399 int target_width = (int)wxRound((double)source.width / factor);
3400 int target_height = (int)wxRound((double)source.height / factor);
3401
3402 int dest_line_length = dest_stride * BPP / 8;
3403
3404 // On MSW, if using DibSections, each scan line starts on a DWORD boundary.
3405 // The DibSection has been allocated to conform with this requirement.
3406#ifdef __PIX_CACHE_DIBSECTION__
3407 dest_line_length = (((dest_stride * 24) + 31) & ~31) >> 3;
3408#endif
3409
3410 if ((target_height == 0) || (target_width == 0)) return false;
3411
3412 // `volatile` may be needed with regard to `sigsetjmp()` below;
3413 // at least it prevents compiler warning about clobbering the variable.
3414 unsigned char *volatile target_data = ppn;
3415 unsigned char *data = ppn;
3416
3417 if (factor > 1) // downsampling
3418 {
3419 if (scale_type == RENDER_HIDEF) {
3420 // Allocate a working buffer based on scale factor
3421 int blur_factor = wxMax(2, Factor);
3422 int wb_size = (source.width) * (blur_factor * 2) * BPP / 8;
3423 s_data = (unsigned char *)malloc(wb_size); // work buffer
3424 unsigned char *pixel;
3425 int y_offset;
3426
3427 for (int y = dest.y; y < (dest.y + dest.height); y++) {
3428 // Read "blur_factor" lines
3429
3430 wxRect s1;
3431 s1.x = source.x;
3432 s1.y = source.y + (int)(y * factor);
3433 s1.width = source.width;
3434 s1.height = blur_factor;
3435 GetChartBits(s1, s_data, 1);
3436
3437 target_data = data + (y * dest_line_length /*dest_stride * BPP/8*/);
3438
3439 for (int x = 0; x < target_width; x++) {
3440 unsigned int avgRed = 0;
3441 unsigned int avgGreen = 0;
3442 unsigned int avgBlue = 0;
3443 unsigned int pixel_count = 0;
3444 unsigned char *pix0 = s_data + BPP / 8 * ((int)(x * factor));
3445 y_offset = 0;
3446
3447 if ((x * Factor) < (Size_X - source.x)) {
3448 // determine average
3449 for (int y1 = 0; y1 < blur_factor; ++y1) {
3450 pixel = pix0 + (BPP / 8 * y_offset);
3451 for (int x1 = 0; x1 < blur_factor; ++x1) {
3452 avgRed += pixel[0];
3453 avgGreen += pixel[1];
3454 avgBlue += pixel[2];
3455
3456 pixel += BPP / 8;
3457
3458 pixel_count++;
3459 }
3460 y_offset += source.width;
3461 }
3462
3463 if (0 == pixel_count) // Protect
3464 pixel_count = 1;
3465
3466 target_data[0] = avgRed / pixel_count; // >> scounter;
3467 target_data[1] = avgGreen / pixel_count; // >> scounter;
3468 target_data[2] = avgBlue / pixel_count; // >> scounter;
3469 target_data += BPP / 8;
3470 } else {
3471 target_data[0] = 0;
3472 target_data[1] = 0;
3473 target_data[2] = 0;
3474 target_data += BPP / 8;
3475 }
3476
3477 } // for x
3478
3479 } // for y
3480
3481 } // SCALE_BILINEAR
3482
3483 else if (scale_type == RENDER_LODEF) {
3484 int get_bits_submap = 1;
3485
3486 int scaler = 16;
3487
3488 if (source.width > 32767) // High underscale can exceed signed math bits
3489 scaler = 8;
3490
3491 int wb_size = (Size_X) * ((/*Factor +*/ 1) * 2) * BPP / 8;
3492 s_data = (unsigned char *)malloc(wb_size); // work buffer
3493
3494 long x_delta = (source.width << scaler) / target_width;
3495 long y_delta = (source.height << scaler) / target_height;
3496
3497 int y = dest.y; // starting here
3498 long ys = dest.y * y_delta;
3499
3500 while (y < dest.y + dest.height) {
3501 // Read 1 line at the right place from the source
3502
3503 wxRect s1;
3504 s1.x = 0;
3505 s1.y = source.y + (ys >> scaler);
3506 s1.width = Size_X;
3507 s1.height = 1;
3508 GetChartBits(s1, s_data, get_bits_submap);
3509
3510 target_data = data + (y * dest_line_length /*dest_stride * BPP/8*/) +
3511 (dest.x * BPP / 8);
3512
3513 long x = (source.x << scaler) + (dest.x * x_delta);
3514 long sizex16 = Size_X << scaler;
3515 int xt = dest.x;
3516
3517 while ((xt < dest.x + dest.width) && (x < 0)) {
3518 target_data[0] = 0;
3519 target_data[1] = 0;
3520 target_data[2] = 0;
3521
3522 target_data += BPP / 8;
3523 x += x_delta;
3524 xt++;
3525 }
3526
3527 while ((xt < dest.x + dest.width) && (x < sizex16)) {
3528 unsigned char *src_pixel = &s_data[(x >> scaler) * BPP / 8];
3529
3530 target_data[0] = src_pixel[0];
3531 target_data[1] = src_pixel[1];
3532 target_data[2] = src_pixel[2];
3533
3534 target_data += BPP / 8;
3535 x += x_delta;
3536 xt++;
3537 }
3538
3539 while (xt < dest.x + dest.width) {
3540 target_data[0] = 0;
3541 target_data[1] = 0;
3542 target_data[2] = 0;
3543
3544 target_data += BPP / 8;
3545 xt++;
3546 }
3547
3548 y++;
3549 ys += y_delta;
3550 }
3551
3552 } // SCALE_SUBSAMP
3553
3554 } else // factor < 1, overzoom
3555 {
3556 int i = 0;
3557 int j = 0;
3558 unsigned char *target_line_start = NULL;
3559 unsigned char *target_data_x = NULL;
3560 int y_offset = 0;
3561
3562#ifdef __WXGTK__
3563 sigaction(SIGSEGV, NULL,
3564 &sa_all_previous); // save existing action for this signal
3565
3566 struct sigaction temp;
3567 sigaction(SIGSEGV, NULL, &temp); // inspect existing action for this signal
3568
3569 temp.sa_handler = catch_signals_chart; // point to my handler
3570 sigemptyset(&temp.sa_mask); // make the blocking set
3571 // empty, so that all
3572 // other signals will be
3573 // unblocked during my handler
3574 temp.sa_flags = 0;
3575 sigaction(SIGSEGV, &temp, NULL);
3576
3577 if (sigsetjmp(env_chart,
3578 1)) // Something in the below code block faulted....
3579 {
3580 sigaction(SIGSEGV, &sa_all_previous, NULL); // reset signal handler
3581
3582 wxString msg;
3583 msg.Printf(" Caught SIGSEGV on GetandScaleData, Factor < 1");
3584 wxLogMessage(msg);
3585
3586 msg.Printf(
3587 " m_raster_scale_factor: %g source.width: %d dest.y: %d "
3588 "dest.x: %d dest.width: %d dest.height: %d ",
3589 m_raster_scale_factor, source.width, dest.y, dest.x, dest.width,
3590 dest.height);
3591 wxLogMessage(msg);
3592
3593 msg.Printf(
3594 " i: %d j: %d dest_stride: %d target_line_start: %p "
3595 "target_data_x: %p y_offset: %d",
3596 i, j, dest_stride, target_line_start, target_data_x, y_offset);
3597 wxLogMessage(msg);
3598
3599 free(s_data);
3600 return true;
3601
3602 }
3603
3604 else
3605#endif
3606 {
3607
3608 double xd, yd;
3609 latlong_to_chartpix(m_vp_render_last.clat, m_vp_render_last.clon, xd, yd);
3610 double xrd =
3611 xd - (m_vp_render_last.pix_width * m_raster_scale_factor / 2);
3612 double yrd =
3613 yd - (m_vp_render_last.pix_height * m_raster_scale_factor / 2);
3614 double x_vernier = (xrd - wxRound(xrd));
3615 double y_vernier = (yrd - wxRound(yrd));
3616 int x_vernier_i = (int)wxRound(x_vernier / m_raster_scale_factor);
3617 int y_vernier_i = (int)wxRound(y_vernier / m_raster_scale_factor);
3618
3619 // Seems safe enough to read all the required data
3620 // Although we must adjust (increase) temporary allocation for negative
3621 // source.x and for vernier
3622 int sx = wxMax(source.x, 0);
3623 s_data = (unsigned char *)malloc((sx + source.width + 2) *
3624 (source.height + 2) * BPP / 8);
3625
3626 wxRect vsource = source;
3627 vsource.height += 2; // get more bits to allow for vernier
3628 vsource.width += 2;
3629 vsource.x -= 1;
3630 vsource.y -= 1;
3631
3632 GetChartBits(vsource, s_data, 1);
3633 unsigned char *source_data = s_data;
3634
3635 j = dest.y;
3636
3637 while (j < dest.y + dest.height) {
3638 y_offset = (int)((j - y_vernier_i) * m_raster_scale_factor) *
3639 vsource.width; // into the source data
3640
3641 target_line_start =
3642 target_data + (j * dest_line_length /*dest_stride * BPP / 8*/);
3643 target_data_x = target_line_start + ((dest.x) * BPP / 8);
3644
3645 i = dest.x;
3646
3647 // Check data bounds to be sure of not overrunning the upstream buffer
3648 if ((target_data_x + (dest.width * BPP / 8)) >
3649 (target_data + data_size)) {
3650 j = dest.y + dest.height;
3651 } else {
3652 while (i < dest.x + dest.width) {
3653 memcpy(target_data_x,
3654 source_data + BPP / 8 *
3655 (y_offset + (int)((i + x_vernier_i) *
3656 m_raster_scale_factor)),
3657 BPP / 8);
3658
3659 target_data_x += BPP / 8;
3660
3661 i++;
3662 }
3663 }
3664
3665 j++;
3666 }
3667 }
3668#ifdef __WXGTK__
3669 sigaction(SIGSEGV, &sa_all_previous, NULL); // reset signal handler
3670#endif
3671 }
3672
3673 free(s_data);
3674
3675 return true;
3676}
3677
3678bool ChartBaseBSB::GetChartBits(wxRect &source, unsigned char *pPix,
3679 int sub_samp) {
3680 wxCriticalSectionLocker locker(m_critSect);
3681
3682 int iy;
3683#define FILL_BYTE 0
3684
3685 // Decode the KAP file RLL stream into image pPix
3686
3687 unsigned char *pCP;
3688 pCP = pPix;
3689
3690 iy = source.y;
3691
3692 while (iy < source.y + source.height) {
3693 if ((iy >= 0) && (iy < Size_Y)) {
3694 if (source.x >= 0) {
3695 if ((source.x + source.width) > Size_X) {
3696 if ((Size_X - source.x) < 0)
3697 memset(pCP, FILL_BYTE, source.width * BPP / 8);
3698 else {
3699 BSBGetScanline(pCP, iy, source.x, Size_X, sub_samp);
3700 memset(pCP + (Size_X - source.x) * BPP / 8, FILL_BYTE,
3701 (source.x + source.width - Size_X) * BPP / 8);
3702 }
3703 } else
3704 BSBGetScanline(pCP, iy, source.x, source.x + source.width, sub_samp);
3705 } else {
3706 if ((source.width + source.x) >= 0) {
3707 // Special case, black on left side
3708 // must ensure that (black fill length % sub_samp) == 0
3709
3710 int xfill_corrected = -source.x + (source.x % sub_samp); //+ve
3711 memset(pCP, FILL_BYTE, (xfill_corrected * BPP / 8));
3712 BSBGetScanline(pCP + (xfill_corrected * BPP / 8), iy, 0,
3713 source.width + source.x, sub_samp);
3714
3715 } else {
3716 memset(pCP, FILL_BYTE, source.width * BPP / 8);
3717 }
3718 }
3719 }
3720
3721 else // requested y is off chart
3722 {
3723 memset(pCP, FILL_BYTE, source.width * BPP / 8);
3724 }
3725
3726 pCP += source.width * BPP / 8 * sub_samp;
3727
3728 iy += sub_samp;
3729 } // while iy
3730
3731 return true;
3732}
3733
3734//-----------------------------------------------------------------------------------------------
3735// BSB File Read Support
3736//-----------------------------------------------------------------------------------------------
3737
3738//-----------------------------------------------------------------------------------------------
3739// ReadBSBHdrLine
3740//
3741// Read and return count of a line of BSB header file
3742//-----------------------------------------------------------------------------------------------
3743
3744int ChartBaseBSB::ReadBSBHdrLine(wxInputStream *ifs, char *buf, int buf_len_max)
3745
3746{
3747 char read_char;
3748 char cr_test;
3749 int line_length = 0;
3750 char *lbuf;
3751
3752 lbuf = buf;
3753
3754 while (!ifs->Eof() && line_length < buf_len_max) {
3755 int c = ifs->GetC();
3756 if (c < 0) break;
3757 read_char = c;
3758 if (0x1A == read_char) {
3759 ifs->Ungetch(read_char);
3760 return (0);
3761 }
3762
3763 if (0 == read_char) // embedded erroneous unicode character?
3764 read_char = 0x20;
3765
3766 // Manage continued lines
3767 if (read_char == 10 || read_char == 13) {
3768 // Check to see if there is an extra CR
3769 cr_test = ifs->GetC();
3770 if (cr_test == 13) cr_test = ifs->GetC(); // skip any extra CR
3771
3772 if (cr_test != 10 && cr_test != 13) ifs->Ungetch(cr_test);
3773 read_char = '\n';
3774 }
3775
3776 // Look for continued lines, indicated by ' ' in first position
3777 if (read_char == '\n') {
3778 cr_test = 0;
3779 cr_test = ifs->GetC();
3780
3781 if (cr_test != ' ') {
3782 ifs->Ungetch(cr_test);
3783 *lbuf = '\0';
3784 return line_length;
3785 }
3786
3787 // Merge out leading spaces
3788 while (cr_test == ' ') cr_test = ifs->GetC();
3789 ifs->Ungetch(cr_test);
3790
3791 // Add a comma
3792 *lbuf = ',';
3793 lbuf++;
3794 }
3795
3796 else {
3797 *lbuf = read_char;
3798 lbuf++;
3799 line_length++;
3800 }
3801
3802 } // while
3803
3804 // Terminate line
3805 if (line_length) *(lbuf - 1) = '\0';
3806
3807 return line_length;
3808}
3809
3810//-----------------------------------------------------------------------
3811// Scan a BSB Scan Line from raw data
3812// Leaving stream pointer at start of next line
3813//-----------------------------------------------------------------------
3814int ChartBaseBSB::BSBScanScanline(wxInputStream *pinStream) {
3815 int nLineMarker, nValueShift, iPixel = 0;
3816 unsigned char byValueMask, byCountMask;
3817 unsigned char byNext;
3818 int coffset;
3819
3820 // if(1)
3821 {
3822 // Read the line number.
3823 nLineMarker = 0;
3824 do {
3825 byNext = pinStream->GetC();
3826 nLineMarker = nLineMarker * 128 + (byNext & 0x7f);
3827 } while ((byNext & 0x80) != 0);
3828
3829 // Setup masking values.
3830 nValueShift = 7 - nColorSize;
3831 byValueMask = (((1 << nColorSize)) - 1) << nValueShift;
3832 byCountMask = (1 << (7 - nColorSize)) - 1;
3833
3834 // Read and simulate expansion of runs.
3835
3836 while (((byNext = pinStream->GetC()) != 0) && (iPixel < Size_X)) {
3837 // int nPixValue;
3838 int nRunCount;
3839 // nPixValue = (byNext & byValueMask) >> nValueShift;
3840
3841 nRunCount = byNext & byCountMask;
3842
3843 while ((byNext & 0x80) != 0) {
3844 byNext = pinStream->GetC();
3845 nRunCount = nRunCount * 128 + (byNext & 0x7f);
3846 }
3847
3848 if (iPixel + nRunCount + 1 > Size_X) nRunCount = Size_X - iPixel - 1;
3849
3850 // Store nPixValue in the destination
3851 // memset(pCL, nPixValue, nRunCount+1);
3852 // pCL += nRunCount+1;
3853 iPixel += nRunCount + 1;
3854 }
3855 coffset = pinStream->TellI();
3856 }
3857
3858 return nLineMarker;
3859}
3860// MSVC compiler makes a bad decision about when to inline (or not) some
3861// intrinsics, like memset(). So,... Here is a little hand-crafted memset()
3862// substitue for known short strings. It will be inlined by MSVC compiler
3863// using /02 settings
3864
3865inline void memset_short(unsigned char *dst, unsigned char cbyte, int count) {
3866#if 0 // def __MSVC__
3867 __asm {
3868 pushf // save Direction flag
3869 cld // set direction "up"
3870 mov edi, dst
3871 mov ecx, count
3872 mov al, cbyte
3873 rep stosb
3874 popf
3875 }
3876#else
3877 memset(dst, cbyte, count);
3878#endif
3879}
3880// could use a larger value for slightly less ram but slower random access,
3881// this is chosen as it is also the opengl tile size so should work well
3882#define TILE_SIZE 512
3883
3884// #define USE_OLD_CACHE // removed this (and simplify code below) once the new
3885// method is verified #define PRINT_TIMINGS // enable for profiling
3886
3887#ifdef PRINT_TIMINGS
3888class OCPNStopWatch {
3889public:
3890 OCPNStopWatch() { Reset(); }
3891 void Reset() { clock_gettime(CLOCK_REALTIME, &tp); }
3892
3893 double Time() {
3894 timespec tp_end;
3895 clock_gettime(CLOCK_REALTIME, &tp_end);
3896 return (tp_end.tv_sec - tp.tv_sec) * 1.e3 +
3897 (tp_end.tv_nsec - tp.tv_nsec) / 1.e6;
3898 }
3899
3900private:
3901 timespec tp;
3902};
3903#endif
3904
3905#define FAIL \
3906 do { \
3907 free(pt->pTileOffset); \
3908 pt->pTileOffset = NULL; \
3909 free(pt->pPix); \
3910 pt->pPix = NULL; \
3911 pt->bValid = false; \
3912 return 0; \
3913 } while (0)
3914
3915//-----------------------------------------------------------------------
3916// Get a BSB Scan Line Using Cache and scan line index if available
3917//-----------------------------------------------------------------------
3918int ChartBaseBSB::BSBGetScanline(unsigned char *pLineBuf, int y, int xs, int xl,
3919 int sub_samp)
3920
3921{
3922 unsigned char *prgb = pLineBuf;
3923 int nValueShift;
3924 unsigned char byValueMask, byCountMask;
3925 unsigned char byNext;
3926 CachedLine *pt = NULL, cached_line;
3927 unsigned char *pCL;
3928 int rgbval;
3929 unsigned char *lp;
3930 int ix = xs;
3931 int pos = 0;
3932
3933 if (bUseLineCache && pLineCache) {
3934 // Is the requested line in the cache, and valid?
3935 pt = &pLineCache[y];
3936 } else {
3937 pt = &cached_line;
3938 pt->bValid = false;
3939 }
3940
3941#ifdef PRINT_TIMINGS
3942 OCPNStopWatch sw;
3943 static double ttime;
3944 static int cnt;
3945 cnt++;
3946#endif
3947
3948 if (!pt->bValid) // not valid, allocate
3949 {
3950 pt->size = pline_table[y + 1] - pline_table[y];
3951
3952#ifdef USE_OLD_CACHE
3953 pt->pPix = (unsigned char *)malloc(Size_X);
3954#else
3955 pt->pTileOffset = (TileOffsetCache *)calloc(
3956 sizeof(TileOffsetCache) * (Size_X / TILE_SIZE + 1), 1);
3957 pt->pPix = (unsigned char *)malloc(pt->size);
3958#endif
3959 if (pline_table[y] == 0 || pline_table[y + 1] == 0) FAIL;
3960
3961 // as of 2015, in wxWidgets buffered streams don't test for a zero seek
3962 // so we check here to possibly avoid this seek with a measured performance
3963 // gain
3964 if (ifs_bitmap->TellI() != pline_table[y] &&
3965 wxInvalidOffset == ifs_bitmap->SeekI(pline_table[y], wxFromStart))
3966 FAIL;
3967
3968#ifdef USE_OLD_CACHE
3969 if (pt->size > ifs_bufsize) {
3970 unsigned char *tmp = ifs_buf;
3971 if (!(ifs_buf = (unsigned char *)realloc(ifs_buf, pt->size))) {
3972 free(tmp);
3973 FAIL;
3974 }
3975 ifs_bufsize = pt->size;
3976 }
3977
3978 lp = ifs_buf;
3979#else
3980 lp = pt->pPix;
3981#endif
3982 ifs_bitmap->Read(lp, pt->size);
3983
3984#ifdef USE_OLD_CACHE
3985 pCL = pt->pPix;
3986#else
3987 if (!bUseLineCache) {
3988 ix = 0;
3989 // skip the line number.
3990 do {
3991 if (pos < pt->size) {
3992 byNext = *lp++;
3993 pos++;
3994 } else {
3995 byNext = 0;
3996 }
3997 } while ((byNext & 0x80) != 0);
3998 goto nocachestart;
3999 }
4000 pCL = ifs_buf;
4001
4002 if (Size_X > ifs_bufsize) {
4003 unsigned char *tmp = ifs_buf;
4004 if (!(ifs_buf = (unsigned char *)realloc(ifs_buf, Size_X))) {
4005 free(tmp);
4006 FAIL;
4007 }
4008 ifs_bufsize = Size_X;
4009 }
4010#endif
4011 // At this point, the unexpanded, raw line is at *lp, and the expansion
4012 // destination is pCL
4013
4014 // skip the line number.
4015 do {
4016 if (pos < pt->size) {
4017 byNext = *lp++;
4018 pos++;
4019 } else {
4020 byNext = 0;
4021 }
4022 } while ((byNext & 0x80) != 0);
4023
4024 // Setup masking values.
4025 nValueShift = 7 - nColorSize;
4026 byValueMask = (((1 << nColorSize)) - 1) << nValueShift;
4027 byCountMask = (1 << (7 - nColorSize)) - 1;
4028
4029 // Read and expand runs.
4030 unsigned int iPixel = 0;
4031
4032#ifndef USE_OLD_CACHE
4033 pt->pTileOffset[0].offset = lp - pt->pPix;
4034 pt->pTileOffset[0].pixel = 0;
4035 unsigned int tileindex = 1, nextTile = TILE_SIZE;
4036#endif
4037 unsigned int nRunCount;
4038 unsigned char *end = pt->pPix + pt->size;
4039 while (iPixel < (unsigned int)Size_X)
4040#ifdef USE_OLD_CACHE
4041 {
4042 nPixValue = (byNext & byValueMask) >> nValueShift;
4043
4044 nRunCount = byNext & byCountMask;
4045
4046 while ((byNext & 0x80) != 0) {
4047 byNext = *lp++;
4048 nRunCount = nRunCount * 128 + (byNext & 0x7f);
4049 }
4050
4051 nRunCount++;
4052
4053 if (iPixel + nRunCount >
4054 (unsigned int)Size_X) // protection against corrupt data
4055 nRunCount = nRunCount - iPixel;
4056
4057 // Store nPixValue in the destination
4058 memset_short(pCL + iPixel, nPixValue, nRunCount);
4059 iPixel += nRunCount;
4060 }
4061#else
4062 // build tile offset table for faster random access
4063 {
4064 if (pos < pt->size) {
4065 byNext = *lp++;
4066 pos++;
4067 } else {
4068 byNext = 0;
4069 }
4070 unsigned char *offset = lp - 1;
4071 if (byNext == 0 || lp == end) {
4072 // finished early...corrupt?
4073 while (tileindex < (unsigned int)Size_X / TILE_SIZE + 1) {
4074 pt->pTileOffset[tileindex].offset = pt->pTileOffset[0].offset;
4075 pt->pTileOffset[tileindex].pixel = 0;
4076 tileindex++;
4077 }
4078 break;
4079 }
4080
4081 nRunCount = byNext & byCountMask;
4082
4083 while ((byNext & 0x80) != 0) {
4084 if (pos < pt->size) {
4085 byNext = *lp++;
4086 pos++;
4087 } else {
4088 byNext = 0;
4089 }
4090 nRunCount = nRunCount * 128 + (byNext & 0x7f);
4091 }
4092
4093 nRunCount++;
4094
4095 if (iPixel + nRunCount >
4096 (unsigned int)Size_X) // protection against corrupt data
4097 nRunCount = Size_X - iPixel;
4098
4099 while (iPixel + nRunCount > nextTile) {
4100 pt->pTileOffset[tileindex].offset = offset - pt->pPix;
4101 pt->pTileOffset[tileindex].pixel = iPixel;
4102 tileindex++;
4103 nextTile += TILE_SIZE;
4104 }
4105 iPixel += nRunCount;
4106 }
4107#endif
4108
4109 pt->bValid = true;
4110 }
4111
4112 // Line is valid, de-reference thru proper pallete directly to target
4113
4114 if (xl > Size_X) xl = Size_X;
4115
4116#ifdef USE_OLD_CACHE
4117 pCL = pt->pPix + xs;
4118
4119 // Optimization for most usual case
4120 if ((BPP == 24) && (1 == sub_samp)) {
4121 ix = xs;
4122 while (ix < xl - 1) {
4123 unsigned char cur_by = *pCL;
4124 rgbval = (int)(pPalette[cur_by]);
4125 while ((ix < xl - 1)) {
4126 if (cur_by != *pCL) break;
4127 *((int *)prgb) = rgbval;
4128 prgb += 3;
4129 pCL++;
4130 ix++;
4131 }
4132 }
4133 } else {
4134 int dest_inc_val_bytes = (BPP / 8) * sub_samp;
4135 ix = xs;
4136 while (ix < xl - 1) {
4137 unsigned char cur_by = *pCL;
4138 rgbval = (int)(pPalette[cur_by]);
4139 while ((ix < xl - 1)) {
4140 if (cur_by != *pCL) break;
4141 *((int *)prgb) = rgbval;
4142 prgb += dest_inc_val_bytes;
4143 pCL += sub_samp;
4144 ix += sub_samp;
4145 }
4146 }
4147 }
4148
4149 // Get the last pixel explicitely
4150 // irrespective of the sub_sampling factor
4151
4152 if (xs < xl - 1) {
4153 unsigned char *pCLast = pt->pPix + (xl - 1);
4154 unsigned char *prgb_last = pLineBuf + ((xl - 1) - xs) * BPP / 8;
4155
4156 rgbval = (int)(pPalette[*pCLast]); // last pixel
4157 unsigned char a = rgbval & 0xff;
4158 *prgb_last++ = a;
4159 a = (rgbval >> 8) & 0xff;
4160 *prgb_last++ = a;
4161 a = (rgbval >> 16) & 0xff;
4162 *prgb_last = a;
4163 }
4164#else
4165 {
4166 int tileindex = xs / TILE_SIZE;
4167 int tileoffset = pt->pTileOffset[tileindex].offset;
4168
4169 lp = pt->pPix + tileoffset;
4170 pos = tileoffset;
4171 ix = pt->pTileOffset[tileindex].pixel;
4172 }
4173
4174nocachestart:
4175 nValueShift = 7 - nColorSize;
4176 byValueMask = (((1 << nColorSize)) - 1) << nValueShift;
4177 byCountMask = (1 << (7 - nColorSize)) - 1;
4178 int nPixValue = 0; // satisfy stupid compiler warning
4179 bool bLastPixValueValid = false;
4180 while (ix < xl - 1) {
4181 if (pos < pt->size) {
4182 byNext = *lp++;
4183 pos++;
4184 } else {
4185 break;
4186 }
4187
4188 nPixValue = (byNext & byValueMask) >> nValueShift;
4189 unsigned int nRunCount;
4190
4191 if (byNext == 0)
4192 nRunCount = xl - ix; // corrupted chart, just run to the end
4193 else {
4194 nRunCount = byNext & byCountMask;
4195 while ((byNext & 0x80) != 0) {
4196 if (pos < pt->size) {
4197 byNext = *lp++;
4198 pos++;
4199 } else {
4200 nRunCount = xl - ix; // corrupted chart, just run to the end
4201 break;
4202 }
4203 nRunCount = nRunCount * 128 + (byNext & 0x7f);
4204 }
4205
4206 nRunCount++;
4207 }
4208
4209 if (ix < xs) {
4210 if (ix + nRunCount <= (unsigned int)xs) {
4211 ix += nRunCount;
4212 continue;
4213 }
4214 nRunCount -= xs - ix;
4215 ix = xs;
4216 }
4217
4218 if (ix + nRunCount >= (unsigned int)xl) {
4219 nRunCount = xl - 1 - ix;
4220 bLastPixValueValid = true;
4221 }
4222
4223 rgbval = (int)(pPalette[nPixValue]);
4224
4225 // Optimization for most usual case
4226 // currently this is the only case possible...
4227 // if((BPP == 24) && (1 == sub_samp))
4228 {
4229 int count = nRunCount;
4230 if (count < 16) {
4231 // for short runs, use simple loop
4232 while (count--) {
4233 *(uint32_t *)prgb = rgbval;
4234 prgb += 3;
4235 }
4236 } else if (rgbval == 0 || rgbval == 0xffffff) {
4237 // optimization for black or white (could work for any gray too)
4238 memset(prgb, rgbval, nRunCount * 3);
4239 prgb += nRunCount * 3;
4240 } else {
4241 // note: this may not be optimal for all processors and compilers
4242 // I optimized for x86_64 using gcc with -O3
4243 // it is probably possible to gain even faster performance by ensuring
4244 // alignment to 16 or 32byte boundary (depending on processor) then
4245 // using inline assembly
4246
4247#ifdef __ARM_ARCH
4248 // ARM needs 8 byte alignment for *(uint64_T *x) = *(uint64_T *y)
4249 // because the compiler will (probably) use the ldrd/strd instuction
4250 // pair. So, advance the prgb pointer until it is 8-byte aligned, and
4251 // then carry on if enough bytes are left to process as 64 bit elements
4252
4253 if ((long)prgb & 7) {
4254 while (count--) {
4255 *(uint32_t *)prgb = rgbval;
4256 prgb += 3;
4257 if (!((long)prgb & 7)) {
4258 if (count >= 8) break;
4259 }
4260 }
4261 }
4262#endif
4263
4264 // fill first 24 bytes
4265 uint64_t *b = (uint64_t *)prgb;
4266 for (int i = 0; i < 8; i++) {
4267 *(uint32_t *)prgb = rgbval;
4268 prgb += 3;
4269 }
4270 count -= 8;
4271
4272 // fill in blocks of 24 bytes
4273 uint64_t *y = (uint64_t *)prgb;
4274 int count_d8 = count >> 3;
4275 prgb += 24 * count_d8;
4276 while (count_d8--) {
4277 *y++ = b[0];
4278 *y++ = b[1];
4279 *y++ = b[2];
4280 }
4281
4282 // fill remaining bytes
4283 int rcount = count & 0x7;
4284 while (rcount--) {
4285 *(uint32_t *)prgb = rgbval;
4286 prgb += 3;
4287 }
4288 }
4289 }
4290
4291 ix += nRunCount;
4292 }
4293
4294 // Get the last pixel explicitely
4295 // irrespective of the sub_sampling factor
4296
4297 if (ix < xl) {
4298 if (!bLastPixValueValid) {
4299 if (pos < pt->size) {
4300 byNext = *lp++;
4301 pos++;
4302 } else {
4303 byNext = 0;
4304 }
4305 nPixValue = (byNext & byValueMask) >> nValueShift;
4306 }
4307 rgbval = (int)(pPalette[nPixValue]); // last pixel
4308 unsigned char a = rgbval & 0xff;
4309
4310 *prgb++ = a;
4311 a = (rgbval >> 8) & 0xff;
4312 *prgb++ = a;
4313 a = (rgbval >> 16) & 0xff;
4314 *prgb = a;
4315 }
4316#endif
4317
4318#ifdef PRINT_TIMINGS
4319 ttime += sw.Time();
4320
4321 if (cnt == 500000) {
4322 static int d;
4323 printf("cache time: %d %f\n", d, ttime / 1000.0);
4324 cnt = 0;
4325 d++;
4326 // ttime = 0;
4327 }
4328#endif
4329
4330 if (!bUseLineCache) {
4331#ifndef USE_OLD_CACHE
4332 free(pt->pTileOffset);
4333#endif
4334 free(pt->pPix);
4335 }
4336
4337 return 1;
4338}
4339
4340int *ChartBaseBSB::GetPalettePtr(BSB_Color_Capability color_index) {
4341 if (pPalettes[color_index]) {
4342 if (palette_direction == PaletteFwd)
4343 return (int *)(pPalettes[color_index]->FwdPalette);
4344 else
4345 return (int *)(pPalettes[color_index]->RevPalette);
4346 } else
4347 return NULL;
4348}
4349
4350PaletteDir ChartBaseBSB::GetPaletteDir() {
4351 // make a pixel cache
4352 PixelCache *pc = new PixelCache(4, 4, BPP);
4353 RGBO r = pc->GetRGBO();
4354 delete pc;
4355
4356 if (r == RGB)
4357 return PaletteFwd;
4358 else
4359 return PaletteRev;
4360}
4361
4362bool ChartBaseBSB::AnalyzeSkew() {
4363 double lonmin = 1000;
4364 double lonmax = -1000;
4365 double latmin = 90.;
4366 double latmax = -90.;
4367
4368 int nlonmin, nlonmax;
4369 nlonmin = 0;
4370 nlonmax = 0;
4371
4372 if (0 == nRefpoint) // bad chart georef...
4373 return (1);
4374
4375 for (int n = 0; n < nRefpoint; n++) {
4376 // Longitude
4377 if (pRefTable[n].lonr > lonmax) {
4378 lonmax = pRefTable[n].lonr;
4379 nlonmax = n;
4380 }
4381 if (pRefTable[n].lonr < lonmin) {
4382 lonmin = pRefTable[n].lonr;
4383 nlonmin = n;
4384 }
4385
4386 // Latitude
4387 if (pRefTable[n].latr < latmin) {
4388 latmin = pRefTable[n].latr;
4389 }
4390 if (pRefTable[n].latr > latmax) {
4391 latmax = pRefTable[n].latr;
4392 }
4393 }
4394
4395 // Special case for charts which cross the IDL
4396 if ((lonmin * lonmax) < 0) {
4397 if (pRefTable[nlonmin].xr > pRefTable[nlonmax].xr) {
4398 // walk the reference table and add 360 to any longitude which is < 0
4399 for (int n = 0; n < nRefpoint; n++) {
4400 if (pRefTable[n].lonr < 0.0) pRefTable[n].lonr += 360.;
4401 }
4402
4403 // And recalculate the min/max
4404 lonmin = 1000;
4405 lonmax = -1000;
4406
4407 for (int n = 0; n < nRefpoint; n++) {
4408 // Longitude
4409 if (pRefTable[n].lonr > lonmax) {
4410 lonmax = pRefTable[n].lonr;
4411 nlonmax = n;
4412 }
4413 if (pRefTable[n].lonr < lonmin) {
4414 lonmin = pRefTable[n].lonr;
4415 nlonmin = n;
4416 }
4417
4418 // Latitude
4419 if (pRefTable[n].latr < latmin) {
4420 latmin = pRefTable[n].latr;
4421 }
4422 if (pRefTable[n].latr > latmax) {
4423 latmax = pRefTable[n].latr;
4424 }
4425 }
4426 }
4427 }
4428
4429 // Find the two REF points that are farthest apart
4430 double dist_max = 0.;
4431 int imax = 0;
4432 int jmax = 0;
4433
4434 for (int i = 0; i < nRefpoint; i++) {
4435 for (int j = i + 1; j < nRefpoint; j++) {
4436 double dx = pRefTable[i].xr - pRefTable[j].xr;
4437 double dy = pRefTable[i].yr - pRefTable[j].yr;
4438 double dist = (dx * dx) + (dy * dy);
4439 if (dist > dist_max) {
4440 dist_max = dist;
4441 imax = i;
4442 jmax = j;
4443 }
4444 }
4445 }
4446
4447 double apparent_skew = 0;
4448
4449 if (m_projection == PROJECTION_MERCATOR) {
4450 double easting0, easting1, northing0, northing1;
4451 // Get the Merc projection of the two REF points
4452 toSM_ECC(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4453 &easting0, &northing0);
4454 toSM_ECC(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4455 &easting1, &northing1);
4456
4457 double skew_proj =
4458 atan2((easting1 - easting0), (northing1 - northing0)) * 180. / PI;
4459 double skew_points = atan2((pRefTable[jmax].yr - pRefTable[imax].yr),
4460 (pRefTable[jmax].xr - pRefTable[imax].xr)) *
4461 180. / PI;
4462
4463 apparent_skew = skew_points - skew_proj + 90.;
4464
4465 // normalize to +/- 180.
4466 if (fabs(apparent_skew) > 180.) {
4467 if (apparent_skew < 0.)
4468 apparent_skew += 360.;
4469 else
4470 apparent_skew -= 360.;
4471 }
4472 }
4473
4474 else if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
4475 double easting0, easting1, northing0, northing1;
4476 // Get the TMerc projection of the two REF points
4477 toTM(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4478 &easting0, &northing0);
4479 toTM(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4480 &easting1, &northing1);
4481
4482 double skew_proj =
4483 atan2((easting1 - easting0), (northing1 - northing0)) * 180. / PI;
4484 double skew_points = atan2((pRefTable[jmax].yr - pRefTable[imax].yr),
4485 (pRefTable[jmax].xr - pRefTable[imax].xr)) *
4486 180. / PI;
4487
4488 apparent_skew = skew_points - skew_proj + 90.;
4489
4490 // normalize to +/- 180.
4491 if (fabs(apparent_skew) > 180.) {
4492 if (apparent_skew < 0.)
4493 apparent_skew += 360.;
4494 else
4495 apparent_skew -= 360.;
4496 }
4497
4498 if (fabs(apparent_skew - m_Chart_Skew) > 2) { // measured skew OK?
4499 // it may be that the projection longitude is simply wrong.
4500 // Check it.
4501 double dskew = fabs(apparent_skew - m_Chart_Skew);
4502 if ((m_proj_lon < lonmin) || (m_proj_lon > lonmax)) {
4503 // Try a projection longitude that is mid-meridian in the chart.
4504 double tentative_proj_lon = (lonmin + lonmax) / 2.;
4505
4506 toTM(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat,
4507 tentative_proj_lon, &easting0, &northing0);
4508 toTM(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat,
4509 tentative_proj_lon, &easting1, &northing1);
4510
4511 skew_proj =
4512 atan2((easting1 - easting0), (northing1 - northing0)) * 180. / PI;
4513 skew_points = atan2((pRefTable[jmax].yr - pRefTable[imax].yr),
4514 (pRefTable[jmax].xr - pRefTable[imax].xr)) *
4515 180. / PI;
4516
4517 apparent_skew = skew_points - skew_proj + 90.;
4518
4519 // normalize to +/- 180.
4520 if (fabs(apparent_skew) > 180.) {
4521 if (apparent_skew < 0.)
4522 apparent_skew += 360.;
4523 else
4524 apparent_skew -= 360.;
4525 }
4526
4527 // Better? If so, adopt the adjusted projection longitude
4528 if (fabs(apparent_skew - m_Chart_Skew) < dskew) {
4529 m_proj_lon = tentative_proj_lon;
4530 }
4531 }
4532 }
4533 } else // For all other projections, assume that skew specified in header is
4534 // correct
4535 apparent_skew = m_Chart_Skew;
4536
4537 if (fabs(apparent_skew - m_Chart_Skew) >
4538 2) { // measured skew is more than 2 degrees
4539 m_Chart_Skew = apparent_skew; // different from stated skew
4540
4541 wxString msg = " Warning: Skew override on chart ";
4542 msg.Append(m_FullPath);
4543 wxString msg1;
4544 msg1.Printf(" is %5g degrees", apparent_skew);
4545 msg.Append(msg1);
4546
4547 wxLogMessage(msg);
4548
4549 return false;
4550 }
4551
4552 return true;
4553}
4554
4555int ChartBaseBSB::AnalyzeRefpoints(bool b_testSolution) {
4556 int i, n;
4557 double elt, elg;
4558
4559 // Calculate the max/min reference points
4560
4561 float lonmin = 1000;
4562 float lonmax = -1000;
4563 float latmin = 90.;
4564 float latmax = -90.;
4565
4566 int plonmin = 100000;
4567 int plonmax = 0;
4568 int platmin = 100000;
4569 int platmax = 0;
4570 int nlonmin = 0, nlonmax = 0;
4571
4572 if (0 == nRefpoint) // bad chart georef...
4573 return (1);
4574
4575 for (n = 0; n < nRefpoint; n++) {
4576 // Longitude
4577 if (pRefTable[n].lonr > lonmax) {
4578 lonmax = pRefTable[n].lonr;
4579 plonmax = (int)pRefTable[n].xr;
4580 nlonmax = n;
4581 }
4582 if (pRefTable[n].lonr < lonmin) {
4583 lonmin = pRefTable[n].lonr;
4584 plonmin = (int)pRefTable[n].xr;
4585 nlonmin = n;
4586 }
4587
4588 // Latitude
4589 if (pRefTable[n].latr < latmin) {
4590 latmin = pRefTable[n].latr;
4591 platmin = (int)pRefTable[n].yr;
4592 }
4593 if (pRefTable[n].latr > latmax) {
4594 latmax = pRefTable[n].latr;
4595 platmax = (int)pRefTable[n].yr;
4596 }
4597 }
4598
4599 // Special case for charts which cross the IDL
4600 if ((lonmin * lonmax) < 0) {
4601 if (pRefTable[nlonmin].xr > pRefTable[nlonmax].xr) {
4602 // walk the reference table and add 360 to any longitude which is < 0
4603 for (n = 0; n < nRefpoint; n++) {
4604 if (pRefTable[n].lonr < 0.0) pRefTable[n].lonr += 360.;
4605 }
4606
4607 // And recalculate the min/max
4608 lonmin = 1000;
4609 lonmax = -1000;
4610
4611 for (n = 0; n < nRefpoint; n++) {
4612 // Longitude
4613 if (pRefTable[n].lonr > lonmax) {
4614 lonmax = pRefTable[n].lonr;
4615 plonmax = (int)pRefTable[n].xr;
4616 nlonmax = n;
4617 }
4618 if (pRefTable[n].lonr < lonmin) {
4619 lonmin = pRefTable[n].lonr;
4620 plonmin = (int)pRefTable[n].xr;
4621 nlonmin = n;
4622 }
4623
4624 // Latitude
4625 if (pRefTable[n].latr < latmin) {
4626 latmin = pRefTable[n].latr;
4627 platmin = (int)pRefTable[n].yr;
4628 }
4629 if (pRefTable[n].latr > latmax) {
4630 latmax = pRefTable[n].latr;
4631 platmax = (int)pRefTable[n].yr;
4632 }
4633 }
4634 }
4635 }
4636
4637 // Build the Control Point Structure, etc
4638 cPoints.count = nRefpoint;
4639 if (cPoints.status) {
4640 // AnalyzeRefpoints can be called twice
4641 free(cPoints.tx);
4642 free(cPoints.ty);
4643 free(cPoints.lon);
4644 free(cPoints.lat);
4645
4646 free(cPoints.pwx);
4647 free(cPoints.wpx);
4648 free(cPoints.pwy);
4649 free(cPoints.wpy);
4650 }
4651
4652 cPoints.tx = (double *)malloc(nRefpoint * sizeof(double));
4653 cPoints.ty = (double *)malloc(nRefpoint * sizeof(double));
4654 cPoints.lon = (double *)malloc(nRefpoint * sizeof(double));
4655 cPoints.lat = (double *)malloc(nRefpoint * sizeof(double));
4656
4657 cPoints.pwx = (double *)malloc(12 * sizeof(double));
4658 cPoints.wpx = (double *)malloc(12 * sizeof(double));
4659 cPoints.pwy = (double *)malloc(12 * sizeof(double));
4660 cPoints.wpy = (double *)malloc(12 * sizeof(double));
4661 cPoints.status = 1;
4662
4663 // Find the two REF points that are farthest apart
4664 double dist_max = 0.;
4665 int imax = 0;
4666 int jmax = 0;
4667
4668 for (i = 0; i < nRefpoint; i++) {
4669 for (int j = i + 1; j < nRefpoint; j++) {
4670 double dx = pRefTable[i].xr - pRefTable[j].xr;
4671 double dy = pRefTable[i].yr - pRefTable[j].yr;
4672 double dist = (dx * dx) + (dy * dy);
4673 if (dist > dist_max) {
4674 dist_max = dist;
4675 imax = i;
4676 jmax = j;
4677 }
4678 }
4679 }
4680
4681 // Georef solution depends on projection type
4682
4683 if (m_projection == PROJECTION_TRANSVERSE_MERCATOR) {
4684 double easting0, easting1, northing0, northing1;
4685 // Get the TMerc projection of the two REF points
4686 toTM(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4687 &easting0, &northing0);
4688 toTM(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4689 &easting1, &northing1);
4690
4691 // Calculate the scale factor using exact REF point math
4692 double dx2 = (pRefTable[jmax].xr - pRefTable[imax].xr) *
4693 (pRefTable[jmax].xr - pRefTable[imax].xr);
4694 double dy2 = (pRefTable[jmax].yr - pRefTable[imax].yr) *
4695 (pRefTable[jmax].yr - pRefTable[imax].yr);
4696 double dn2 = (northing1 - northing0) * (northing1 - northing0);
4697 double de2 = (easting1 - easting0) * (easting1 - easting0);
4698
4699 m_ppm_avg = sqrt(dx2 + dy2) / sqrt(dn2 + de2);
4700
4701 // Set up and solve polynomial solution for pix<->east/north as projected
4702 // Fill the cpoints structure with pixel points and transformed lat/lon
4703
4704 for (int n = 0; n < nRefpoint; n++) {
4705 double easting, northing;
4706 toTM(pRefTable[n].latr, pRefTable[n].lonr, m_proj_lat, m_proj_lon,
4707 &easting, &northing);
4708
4709 cPoints.tx[n] = pRefTable[n].xr;
4710 cPoints.ty[n] = pRefTable[n].yr;
4711 cPoints.lon[n] = easting;
4712 cPoints.lat[n] = northing;
4713 }
4714
4715 // Helper parameters
4716 cPoints.txmax = plonmax;
4717 cPoints.txmin = plonmin;
4718 cPoints.tymax = platmax;
4719 cPoints.tymin = platmin;
4720 toTM(latmax, lonmax, m_proj_lat, m_proj_lon, &cPoints.lonmax,
4721 &cPoints.latmax);
4722 toTM(latmin, lonmin, m_proj_lat, m_proj_lon, &cPoints.lonmin,
4723 &cPoints.latmin);
4724
4725 Georef_Calculate_Coefficients_Proj(&cPoints);
4726
4727 }
4728
4729 else if (m_projection == PROJECTION_MERCATOR) {
4730 double easting0, easting1, northing0, northing1;
4731 // Get the Merc projection of the two REF points
4732 toSM_ECC(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4733 &easting0, &northing0);
4734 toSM_ECC(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4735 &easting1, &northing1);
4736
4737 // Calculate the scale factor using exact REF point math
4738 // double dx = (pRefTable[jmax].xr - pRefTable[imax].xr);
4739 // double de = (easting1 - easting0);
4740 // m_ppm_avg = fabs(dx / de);
4741
4742 double dx2 = (pRefTable[jmax].xr - pRefTable[imax].xr) *
4743 (pRefTable[jmax].xr - pRefTable[imax].xr);
4744 double dy2 = (pRefTable[jmax].yr - pRefTable[imax].yr) *
4745 (pRefTable[jmax].yr - pRefTable[imax].yr);
4746 double dn2 = (northing1 - northing0) * (northing1 - northing0);
4747 double de2 = (easting1 - easting0) * (easting1 - easting0);
4748
4749 m_ppm_avg = sqrt(dx2 + dy2) / sqrt(dn2 + de2);
4750
4751 // Set up and solve polynomial solution for pix<->east/north as projected
4752 // Fill the cpoints structure with pixel points and transformed lat/lon
4753
4754 for (int n = 0; n < nRefpoint; n++) {
4755 double easting, northing;
4756 toSM_ECC(pRefTable[n].latr, pRefTable[n].lonr, m_proj_lat, m_proj_lon,
4757 &easting, &northing);
4758
4759 cPoints.tx[n] = pRefTable[n].xr;
4760 cPoints.ty[n] = pRefTable[n].yr;
4761 cPoints.lon[n] = easting;
4762 cPoints.lat[n] = northing;
4763 // printf(" x: %g y: %g east: %g north:
4764 // %g\n",pRefTable[n].xr, pRefTable[n].yr, easting,
4765 // northing);
4766 }
4767
4768 // Helper parameters
4769 cPoints.txmax = plonmax;
4770 cPoints.txmin = plonmin;
4771 cPoints.tymax = platmax;
4772 cPoints.tymin = platmin;
4773 toSM_ECC(latmax, lonmax, m_proj_lat, m_proj_lon, &cPoints.lonmax,
4774 &cPoints.latmax);
4775 toSM_ECC(latmin, lonmin, m_proj_lat, m_proj_lon, &cPoints.lonmin,
4776 &cPoints.latmin);
4777
4778 Georef_Calculate_Coefficients_Proj(&cPoints);
4779
4780 // for(int h=0 ; h < 10 ; h++)
4781 // printf("pix to east %d %g\n", h, cPoints.pwx[h]); //
4782 // pix to lon
4783 // for(int h=0 ; h < 10 ; h++)
4784 // printf("east to pix %d %g\n", h, cPoints.wpx[h]); //
4785 // lon to pix
4786
4787 /*
4788 if ((0 != m_Chart_DU ) && (0 != m_Chart_Scale))
4789 {
4790 double m_ppm_avg1 = m_Chart_DU * 39.37 / m_Chart_Scale;
4791 m_ppm_avg1 *= cos(m_proj_lat * PI / 180.); // correct to
4792 chart centroid
4793
4794 printf("BSB chart ppm_avg:%g ppm_avg1:%g\n", m_ppm_avg,
4795 m_ppm_avg1); m_ppm_avg = m_ppm_avg1;
4796 }
4797 */
4798 }
4799
4800 else if (m_projection == PROJECTION_POLYCONIC) {
4801 // This is interesting
4802 // On some BSB V 1.0 Polyconic charts (e.g. 14500_1, 1995), the projection
4803 // parameter Which is taken to be the central meridian of the projection
4804 // is of the wrong sign....
4805
4806 // We check for this case, and make a correction if necessary.....
4807 // Obviously, the projection meridian should be on the chart, i.e. between
4808 // the min and max longitudes....
4809 double proj_meridian = m_proj_lon;
4810
4811 if ((pRefTable[nlonmax].lonr >= -proj_meridian) &&
4812 (-proj_meridian >= pRefTable[nlonmin].lonr))
4813 m_proj_lon = -m_proj_lon;
4814
4815 double easting0, easting1, northing0, northing1;
4816 // Get the Poly projection of the two REF points
4817 toPOLY(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4818 &easting0, &northing0);
4819 toPOLY(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4820 &easting1, &northing1);
4821
4822 // Calculate the scale factor using exact REF point math
4823 double dx2 = (pRefTable[jmax].xr - pRefTable[imax].xr) *
4824 (pRefTable[jmax].xr - pRefTable[imax].xr);
4825 double dy2 = (pRefTable[jmax].yr - pRefTable[imax].yr) *
4826 (pRefTable[jmax].yr - pRefTable[imax].yr);
4827 double dn2 = (northing1 - northing0) * (northing1 - northing0);
4828 double de2 = (easting1 - easting0) * (easting1 - easting0);
4829
4830 m_ppm_avg = sqrt(dx2 + dy2) / sqrt(dn2 + de2);
4831
4832 // Sanity check
4833 // double ref_dist = DistGreatCircle(pRefTable[imax].latr,
4834 // pRefTable[imax].lonr, pRefTable[jmax].latr,
4835 // pRefTable[jmax].lonr); ref_dist *= 1852; //To Meters double
4836 // ref_dist_transform = sqrt(dn2 + de2); //Also meters
4837 // double error = (ref_dist - ref_dist_transform)/ref_dist;
4838
4839 // Set up and solve polynomial solution for pix<->cartesian east/north as
4840 // projected
4841 // Fill the cpoints structure with pixel points and transformed lat/lon
4842
4843 for (int n = 0; n < nRefpoint; n++) {
4844 double easting, northing;
4845 toPOLY(pRefTable[n].latr, pRefTable[n].lonr, m_proj_lat, m_proj_lon,
4846 &easting, &northing);
4847
4848 // Round trip check for debugging....
4849 // double lat, lon;
4850 // fromPOLY(easting, northing, m_proj_lat, m_proj_lon,
4851 // &lat, &lon);
4852
4853 cPoints.tx[n] = pRefTable[n].xr;
4854 cPoints.ty[n] = pRefTable[n].yr;
4855 cPoints.lon[n] = easting;
4856 cPoints.lat[n] = northing;
4857 // printf(" x: %g y: %g east: %g north:
4858 // %g\n",pRefTable[n].xr, pRefTable[n].yr, easting,
4859 // northing);
4860 }
4861
4862 // Helper parameters
4863 cPoints.txmax = plonmax;
4864 cPoints.txmin = plonmin;
4865 cPoints.tymax = platmax;
4866 cPoints.tymin = platmin;
4867 toPOLY(latmax, lonmax, m_proj_lat, m_proj_lon, &cPoints.lonmax,
4868 &cPoints.latmax);
4869 toPOLY(latmin, lonmin, m_proj_lat, m_proj_lon, &cPoints.lonmin,
4870 &cPoints.latmin);
4871
4872 Georef_Calculate_Coefficients_Proj(&cPoints);
4873
4874 // for(int h=0 ; h < 10 ; h++)
4875 // printf("pix to east %d %g\n", h, cPoints.pwx[h]); //
4876 // pix to lon
4877 // for(int h=0 ; h < 10 ; h++)
4878 // printf("east to pix %d %g\n", h, cPoints.wpx[h]); //
4879 // lon to pix
4880
4881 }
4882
4883 // Any other projection had better have embedded coefficients
4884 else if (bHaveEmbeddedGeoref) {
4885 // Use a Mercator Projection to get a rough ppm.
4886 double easting0, easting1, northing0, northing1;
4887 // Get the Merc projection of the two REF points
4888 toSM_ECC(pRefTable[imax].latr, pRefTable[imax].lonr, m_proj_lat, m_proj_lon,
4889 &easting0, &northing0);
4890 toSM_ECC(pRefTable[jmax].latr, pRefTable[jmax].lonr, m_proj_lat, m_proj_lon,
4891 &easting1, &northing1);
4892
4893 // Calculate the scale factor using exact REF point math in x(longitude)
4894 // direction
4895
4896 double dx = (pRefTable[jmax].xr - pRefTable[imax].xr);
4897 double de = (easting1 - easting0);
4898
4899 m_ppm_avg = fabs(dx / de);
4900
4901 m_ExtraInfo = "---<<< Warning: Chart georef accuracy may be poor. >>>---";
4902 }
4903
4904 else
4905 m_ppm_avg = 1.0; // absolute fallback to prevent div-0 errors
4906
4907#if 0
4908 // Alternate Skew verification
4909 ViewPort vps;
4910 vps.clat = pRefTable[0].latr;
4911 vps.clon = pRefTable[0].lonr;
4912 vps.view_scale_ppm = m_ppm_avg;
4913 vps.skew = 0.;
4914 vps.pix_width = 1000;
4915 vps.pix_height = 1000;
4916
4917 int x1, y1, x2, y2;
4918 latlong_to_pix_vp(latmin, (lonmax + lonmin)/2., x1, y1, vps);
4919 latlong_to_pix_vp(latmax, (lonmax + lonmin)/2., x2, y2, vps);
4920
4921 double apparent_skew = (atan2( (y2-y1), (x2-x1) ) * 180./PI) + 90.;
4922 if(apparent_skew < 0.)
4923 apparent_skew += 360;
4924 if(apparent_skew > 360.)
4925 apparent_skew -= 360;
4926
4927 if(fabs( apparent_skew - m_Chart_Skew ) > 2) { // measured skew is more than 2 degress different
4928 m_Chart_Skew = apparent_skew;
4929 }
4930#endif
4931
4932 if (!b_testSolution) return (0);
4933
4934 // Do a last little test using a synthetic ViewPort of nominal size.....
4935 ViewPort vp;
4936 vp.clat = pRefTable[0].latr;
4937 vp.clon = pRefTable[0].lonr;
4938 vp.view_scale_ppm = m_ppm_avg;
4939 vp.skew = 0.;
4940 vp.pix_width = 1000;
4941 vp.pix_height = 1000;
4942 // vp.rv_rect = wxRect(0,0, vp.pix_width, vp.pix_height);
4943 SetVPRasterParms(vp);
4944
4945 double xpl_err_max = 0;
4946 double ypl_err_max = 0;
4947 int px, py;
4948
4949 int pxref, pyref;
4950 pxref = (int)pRefTable[0].xr;
4951 pyref = (int)pRefTable[0].yr;
4952
4953 for (i = 0; i < nRefpoint; i++) {
4954 px = (int)(vp.pix_width / 2 + pRefTable[i].xr) - pxref;
4955 py = (int)(vp.pix_height / 2 + pRefTable[i].yr) - pyref;
4956
4957 vp_pix_to_latlong(vp, px, py, &elt, &elg);
4958
4959 double lat_error = elt - pRefTable[i].latr;
4960 pRefTable[i].ypl_error = lat_error;
4961
4962 double lon_error = elg - pRefTable[i].lonr;
4963
4964 // Longitude errors could be compounded by prior adjustment to ref points
4965 if (fabs(lon_error) > 180.) {
4966 if (lon_error > 0.)
4967 lon_error -= 360.;
4968 else if (lon_error < 0.)
4969 lon_error += 360.;
4970 }
4971 pRefTable[i].xpl_error = lon_error;
4972
4973 if (fabs(pRefTable[i].ypl_error) > fabs(ypl_err_max))
4974 ypl_err_max = pRefTable[i].ypl_error;
4975 if (fabs(pRefTable[i].xpl_error) > fabs(xpl_err_max))
4976 xpl_err_max = pRefTable[i].xpl_error;
4977 }
4978
4979 Chart_Error_Factor = fmax(fabs(xpl_err_max / (lonmax - lonmin)),
4980 fabs(ypl_err_max / (latmax - latmin)));
4981 double chart_error_meters =
4982 fmax(fabs(xpl_err_max * 60. * 1852.), fabs(ypl_err_max * 60. * 1852.));
4983 // calculate a nominal pixel error
4984 // Assume a modern display has about 4000 pixels/meter.
4985 // Assume the chart is to be displayed at nominal printed scale
4986 double chart_error_pixels = chart_error_meters * 4000. / m_Chart_Scale;
4987
4988 // Good enough for navigation?
4989 int max_pixel_error = 4;
4990
4991 if (chart_error_pixels > max_pixel_error) {
4992 wxString msg =
4993 " VP Final Check: Georeference Chart_Error_Factor on chart ";
4994 msg.Append(m_FullPath);
4995 wxString msg1;
4996 msg1.Printf(" is %5g \n nominal pixel error is: %5g",
4997 Chart_Error_Factor, chart_error_pixels);
4998 msg.Append(msg1);
4999
5000 wxLogMessage(msg);
5001
5002 m_ExtraInfo = "---<<< Warning: Chart georef accuracy is poor. >>>---";
5003 }
5004
5005 // Try again with my calculated georef
5006 // This problem was found on NOAA 514_1.KAP. The embedded coefficients are
5007 // just wrong....
5008 if ((chart_error_pixels > max_pixel_error) && bHaveEmbeddedGeoref) {
5009 wxString msg =
5010 " Trying again with internally calculated georef solution ";
5011 wxLogMessage(msg);
5012
5013 bHaveEmbeddedGeoref = false;
5014 SetVPRasterParms(vp);
5015
5016 xpl_err_max = 0;
5017 ypl_err_max = 0;
5018
5019 pxref = (int)pRefTable[0].xr;
5020 pyref = (int)pRefTable[0].yr;
5021
5022 for (i = 0; i < nRefpoint; i++) {
5023 px = (int)(vp.pix_width / 2 + pRefTable[i].xr) - pxref;
5024 py = (int)(vp.pix_height / 2 + pRefTable[i].yr) - pyref;
5025
5026 vp_pix_to_latlong(vp, px, py, &elt, &elg);
5027
5028 double lat_error = elt - pRefTable[i].latr;
5029 pRefTable[i].ypl_error = lat_error;
5030
5031 double lon_error = elg - pRefTable[i].lonr;
5032
5033 // Longitude errors could be compounded by prior adjustment to ref points
5034 if (fabs(lon_error) > 180.) {
5035 if (lon_error > 0.)
5036 lon_error -= 360.;
5037 else if (lon_error < 0.)
5038 lon_error += 360.;
5039 }
5040 pRefTable[i].xpl_error = lon_error;
5041
5042 if (fabs(pRefTable[i].ypl_error) > fabs(ypl_err_max))
5043 ypl_err_max = pRefTable[i].ypl_error;
5044 if (fabs(pRefTable[i].xpl_error) > fabs(xpl_err_max))
5045 xpl_err_max = pRefTable[i].xpl_error;
5046 }
5047
5048 Chart_Error_Factor = fmax(fabs(xpl_err_max / (lonmax - lonmin)),
5049 fabs(ypl_err_max / (latmax - latmin)));
5050
5051 chart_error_meters =
5052 fmax(fabs(xpl_err_max * 60. * 1852.), fabs(ypl_err_max * 60. * 1852.));
5053 chart_error_pixels = chart_error_meters * 4000. / m_Chart_Scale;
5054
5055 // Good enough for navigation?
5056 if (chart_error_pixels > max_pixel_error) {
5057 wxString msg =
5058 " VP Final Check with internal georef: Georeference "
5059 "Chart_Error_Factor on chart ";
5060 msg.Append(m_FullPath);
5061 wxString msg1;
5062 msg1.Printf(" is %5g\n nominal pixel error is: %5g",
5063 Chart_Error_Factor, chart_error_pixels);
5064 msg.Append(msg1);
5065
5066 wxLogMessage(msg);
5067
5068 m_ExtraInfo = "---<<< Warning: Chart georef accuracy is poor. >>>---";
5069 } else {
5070 wxString msg = " Result: OK, Internal georef solution used.";
5071
5072 wxLogMessage(msg);
5073
5074 m_ExtraInfo = "";
5075 }
5076 }
5077
5078 return (0);
5079}
5080
5081double ChartBaseBSB::AdjustLongitude(double lon) {
5082 double lond = (m_LonMin + m_LonMax) / 2 - lon;
5083 if (lond > 180)
5084 return lon + 360;
5085 else if (lond < -180)
5086 return lon - 360;
5087 return lon;
5088}
5089
5090/*
5091 * Extracted from bsb_io.c - implementation of libbsb reading and writing
5092 *
5093 * Copyright (C) 2000 Stuart Cunningham <stuart_hc@users.sourceforge.net>
5094 *
5095 * This library is free software; you can redistribute it and/or
5096 * modify it under the terms of the GNU Lesser General Public
5097 * License as published by the Free Software Foundation; either
5098 * version 2.1 of the License, or (at your option) any later version.
5099 *
5100 * This library is distributed in the hope that it will be useful,
5101 * but WITHOUT ANY WARRANTY; without even the implied warranty of
5102 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
5103 * Lesser General Public License for more details.
5104 *
5105 * You should have received a copy of the GNU Lesser General Public
5106 * License along with this library; if not, write to the Free Software
5107 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
5108 *
5109 *
5110 */
5111
5121static double polytrans(double *coeff, double lon, double lat) {
5122 double ret = coeff[0] + coeff[1] * lon + coeff[2] * lat;
5123 ret += coeff[3] * lon * lon;
5124 ret += coeff[4] * lon * lat;
5125 ret += coeff[5] * lat * lat;
5126 ret += coeff[6] * lon * lon * lon;
5127 ret += coeff[7] * lon * lon * lat;
5128 ret += coeff[8] * lon * lat * lat;
5129 ret += coeff[9] * lat * lat * lat;
5130 // ret += coeff[10]*lat*lat*lat*lat;
5131 // ret += coeff[11]*lat*lat*lat*lat*lat;
5132
5133 // for(int n = 0 ; n < 10 ; n++)
5134 // printf(" %g\n", coeff[n]);
5135
5136 return ret;
5137}
5138
5139#if 0
5151extern int bsb_LLtoXY(BSBImage *p, double lon, double lat, int* x, int* y)
5152{
5153 /* change longitude phase (CPH) */
5154 lon = (lon < 0) ? lon + p->cph : lon - p->cph;
5155 double xd = polytrans( p->wpx, lon, lat );
5156 double yd = polytrans( p->wpy, lon, lat );
5157 *x = (int)(xd + 0.5);
5158 *y = (int)(yd + 0.5);
5159 return 1;
5160}
5161
5173extern int bsb_XYtoLL(BSBImage *p, int x, int y, double* lonout, double* latout)
5174{
5175 double lon = polytrans( p->pwx, x, y );
5176 lon = (lon < 0) ? lon + p->cph : lon - p->cph;
5177 *lonout = lon;
5178 *latout = polytrans( p->pwy, x, y );
5179 return 1;
5180}
5181
5182#endif
XZ compressed charts support.
BSB chart management.
virtual double GetNearestPreferredScalePPM(double target_scale_ppm)
Find the nearest preferred viewport scale (in pixels/meter) for this chart.
double GetClosestValidNaturalScalePPM(double target_scale, double scale_factor_min, double scale_factor_max)
Find closest valid scale that's a power of 2 multiple of chart's native scale.
An iterator class for OCPNRegion.
A wrapper class for wxRegion with additional functionality.
Definition ocpn_region.h:37
ViewPort - Core geographic projection and coordinate transformation engine.
Definition viewport.h:56
double view_scale_ppm
Requested view scale in physical pixels per meter (ppm), before applying projections.
Definition viewport.h:204
int pix_height
Height of the viewport in physical pixels.
Definition viewport.h:233
int pix_width
Width of the viewport in physical pixels.
Definition viewport.h:231
double skew
Angular distortion (shear transform) applied to the viewport in radians.
Definition viewport.h:212
double clon
Center longitude of the viewport in degrees.
Definition viewport.h:199
double clat
Center latitude of the viewport in degrees.
Definition viewport.h:197
Utility functions.
Optimized wxBitmap Object.
double lonmin
Minimum longitude in reference data.
Definition georef.h:83
int count
Number of reference points used.
Definition georef.h:66
double lonmax
Maximum longitude in reference data.
Definition georef.h:82
int txmin
Minimum x value in target space.
Definition georef.h:80
double * pwx
Polynomial coefficients for pixel-to-world longitude transformation.
Definition georef.h:72
double * tx
Array of x-coordinates in target (typically pixel) space.
Definition georef.h:68
double * wpx
Polynomial coefficients for world-to-pixel x transformation.
Definition georef.h:76
double * lat
Array of latitudes corresponding to reference points.
Definition georef.h:71
double * lon
Array of longitudes corresponding to reference points.
Definition georef.h:70
int status
Status of the georeferencing (0 = valid, other values indicate errors)
Definition georef.h:64
int tymin
Minimum y value in target space.
Definition georef.h:81
double * ty
Array of y-coordinates in target (typically pixel) space.
Definition georef.h:69
int tymax
Maximum y value in target space.
Definition georef.h:79
double latmin
Minimum latitude in reference data.
Definition georef.h:85
int txmax
Maximum x value in target space.
Definition georef.h:78
double latmax
Maximum latitude in reference data.
Definition georef.h:84
double * wpy
Polynomial coefficients for world-to-pixel y transformation.
Definition georef.h:77
double * pwy
Polynomial coefficients for pixel-to-world latitude transformation.
Definition georef.h:74