1 /**********************************************************************
2 zyGrib: meteorological GRIB file viewer
3 Copyright (C) 2008 - Jacques Zaninetti - http://www.zygrib.org
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 ***********************************************************************/
18
19 #include "wx/wxprec.h"
20
21 #ifndef WX_PRECOMP
22 #include "wx/wx.h"
23 #endif //precompiled headers
24
25 //#include "dychart.h" // for some compile time fixups
26 //#include "chcanv.h"
27 //#include "cutil.h"
28 //#include "georef.h"
29 #include <wx/graphics.h>
30
31 #include "IsoLine.h"
32 #include "GribSettingsDialog.h"
33 #include "GribOverlayFactory.h"
34
35 #ifdef __OCPN__ANDROID__
36 #include "qdebug.h"
37 #endif
38
39 //static void GenerateSpline(int n, wxPoint points[]);
40 //static void ClearSplineList();
41 wxList ocpn_wx_spline_point_list;
42
43 #include <wx/listimpl.cpp>
44 WX_DEFINE_LIST(MySegList);
45 WX_DEFINE_LIST(MySegListList);
46
47 #ifndef PI
48 #define PI 3.14159
49 #endif
50 #define CTRUE -1
51 #define CFALSE 0
52
53
54 /* Local variables for cohen_sutherland_line_clip: */
55 struct LOC_cohen_sutherland_line_clip {
56 double xmin, xmax, ymin, ymax;
57 } ;
CompOutCode(double x,double y,outcode * code,struct LOC_cohen_sutherland_line_clip * LINK)58 void CompOutCode (double x, double y, outcode *code, struct LOC_cohen_sutherland_line_clip *LINK)
59 {
60 /*Compute outcode for the point (x,y) */
61 *code = 0;
62 if (y > LINK->ymax)
63 *code = 1L << ((long)TOP);
64 else if (y < LINK->ymin)
65 *code = 1L << ((long)BOTTOM);
66 if (x > LINK->xmax)
67 *code |= 1L << ((long)RIGHT);
68 else if (x < LINK->xmin)
69 *code |= 1L << ((long)LEFT);
70 }
71
72
cohen_sutherland_line_clip_d(double * x0,double * y0,double * x1,double * y1,double xmin_,double xmax_,double ymin_,double ymax_)73 ClipResult cohen_sutherland_line_clip_d (double *x0, double *y0, double *x1, double *y1,
74 double xmin_, double xmax_, double ymin_, double ymax_)
75 {
76 /* Cohen-Sutherland clipping algorithm for line P0=(x1,y0) to P1=(x1,y1)
77 * and clip rectangle with diagonal from (xmin,ymin) to (xmax,ymax).*/
78 struct LOC_cohen_sutherland_line_clip V;
79 int done = CFALSE;
80 ClipResult clip = Visible;
81 outcode outcode0, outcode1, outcodeOut;
82 /*Outcodes for P0,P1, and whichever point lies outside the clip rectangle*/
83 double x=0., y=0.;
84
85 V.xmin = xmin_;
86 V.xmax = xmax_;
87 V.ymin = ymin_;
88 V.ymax = ymax_;
89 CompOutCode(*x0, *y0, &outcode0, &V);
90 CompOutCode(*x1, *y1, &outcode1, &V);
91 do {
92 if (outcode0 == 0 && outcode1 == 0) { /*Trivial accept and exit*/
93 done = CTRUE;
94 } else if ((outcode0 & outcode1) != 0) {
95 clip = Invisible;
96 done = CTRUE;
97 }
98 /*Logical intersection is true, so trivial reject and exit.*/
99 else {
100 clip = Visible;
101 /*Failed both tests, so calculate the line segment to clip;
102 * from an outside point to an intersection with clip edge.*/
103 /*At least one endpoint is outside the clip rectangle; pick it.*/
104 if (outcode0 != 0)
105 outcodeOut = outcode0;
106 else
107 outcodeOut = outcode1;
108 /*Now find intersection point;
109 * use formulas y=y0+slope*(x-x0),x=x0+(1/slope)*(y-y0).*/
110
111 if (((1L << ((long)TOP)) & outcodeOut) != 0) {
112 /*Divide line at top of clip rectangle*/
113 x = *x0 + (*x1 - *x0) * (V.ymax - *y0) / (*y1 - *y0);
114 y = V.ymax;
115 } else if (((1L << ((long)BOTTOM)) & outcodeOut) != 0) {
116 /*Divide line at bottom of clip rectangle*/
117 x = *x0 + (*x1 - *x0) * (V.ymin - *y0) / (*y1 - *y0);
118 y = V.ymin;
119 } else if (((1L << ((long)RIGHT)) & outcodeOut) != 0) {
120 /*Divide line at right edge of clip rectangle*/
121 y = *y0 + (*y1 - *y0) * (V.xmax - *x0) / (*x1 - *x0);
122 x = V.xmax;
123 } else if (((1L << ((long)LEFT)) & outcodeOut) != 0) {
124 /*Divide line at left edge of clip rectangle*/
125 y = *y0 + (*y1 - *y0) * (V.xmin - *x0) / (*x1 - *x0);
126 x = V.xmin;
127 }
128 /*Now we move outside point to intersection point to clip,
129 * and get ready for next pass.*/
130 if (outcodeOut == outcode0) {
131 *x0 = x;
132 *y0 = y;
133 CompOutCode(*x0, *y0, &outcode0, &V);
134 } else {
135 *x1 = x;
136 *y1 = y;
137 CompOutCode(*x1, *y1, &outcode1, &V);
138 }
139 }
140 } while (!done);
141 return clip;
142 }
143
cohen_sutherland_line_clip_i(int * x0_,int * y0_,int * x1_,int * y1_,int xmin_,int xmax_,int ymin_,int ymax_)144 ClipResult cohen_sutherland_line_clip_i (int *x0_, int *y0_, int *x1_, int *y1_,
145 int xmin_, int xmax_, int ymin_, int ymax_)
146 {
147 ClipResult ret;
148 double x0,y0,x1,y1;
149 x0 = *x0_;
150 y0 = *y0_;
151 x1 = *x1_;
152 y1 = *y1_;
153 ret = cohen_sutherland_line_clip_d (&x0, &y0, &x1, &y1,
154 (double)xmin_, (double)xmax_,
155 (double)ymin_, (double)ymax_);
156 *x0_ = (int)x0;
157 *y0_ = (int)y0;
158 *x1_ = (int)x1;
159 *y1_ = (int)y1;
160 return ret;
161 }
162
163
round_msvc(double x)164 double round_msvc (double x)
165 {
166 return(floor(x + 0.5));
167
168 }
169
170 //---------------------------------------------------------------
IsoLine(double val,double coeff,double offset,const GribRecord * rec_)171 IsoLine::IsoLine(double val, double coeff, double offset, const GribRecord *rec_)
172 {
173 if(wxGetDisplaySize().x > 0){
174 m_pixelMM = PlugInGetDisplaySizeMM() / wxGetDisplaySize().x;
175 m_pixelMM = wxMax(.02, m_pixelMM); // protect against bad data
176 }
177 else
178 m_pixelMM = 0.27; // semi-standard number...
179
180 value = val/coeff-offset;
181
182 rec = rec_;
183 W = rec_->getNi();
184 H = rec_->getNj();
185
186 //---------------------------------------------------------
187 // Génère la liste des segments.
188 extractIsoLine(rec_);
189
190 value = val;
191
192 if(trace.size() == 0)
193 return;
194
195 // Join the isoline segments into a nice list
196 // Which is end-to-end continuous and unidirectional
197
198 // Create a master wxList of the trace list
199 std::list<Segment *>::iterator it;
200 for (it=trace.begin(); it!=trace.end(); it++)
201 {
202 Segment *seg = *it;
203 seg->bUsed = false;
204 m_seglist.Append(*it);
205 }
206
207 // Isoline may be discontinuous....
208 // So build a list of continuous segments
209 bool bdone = false;
210 while(!bdone)
211 {
212 MySegList *ps = BuildContinuousSegment();
213
214 m_SegListList.Append(ps);
215
216 MySegList::Node *node;
217 Segment *seg;
218
219 // recreate the master list, removing used segs
220
221 node = m_seglist.GetFirst();
222 while(node)
223 {
224 seg = node->GetData();
225 if(seg->bUsed)
226 {
227 m_seglist.Erase(node);
228 node = m_seglist.GetFirst();
229 }
230 else
231 node = node->GetNext();
232 }
233
234 if(0 == m_seglist.GetCount())
235 bdone = true;
236 }
237
238 ///printf("create Isobar : press=%4.0f long=%d\n", pressure/100, trace.size());
239 }
240 //---------------------------------------------------------------
~IsoLine()241 IsoLine::~IsoLine()
242 {
243 //printf("delete Isobar : press=%4.0f long=%d\n", pressure/100, trace.size());
244
245 std::list<Segment *>::iterator it;
246 for (it=trace.begin(); it!=trace.end(); it++) {
247 delete *it;
248 *it = NULL;
249 }
250 trace.clear();
251
252 m_SegListList.DeleteContents(true);
253 m_SegListList.Clear();
254
255 }
256
257
BuildContinuousSegment(void)258 MySegList *IsoLine::BuildContinuousSegment(void)
259 {
260 MySegList::Node *node;
261 Segment *seg;
262
263 MySegList *ret_list = new MySegList;
264
265 // Build a chain extending from the "2" end of the target segment
266 // The joined list, side 2...
267 MySegList segjoin2;
268
269 // Add any first segment to the list
270 node = m_seglist.GetFirst();
271 Segment *seg0 = node->GetData();
272 seg0->bUsed = true;
273 segjoin2.Append(seg0);
274
275 Segment *tseg = seg0;
276
277 while(tseg)
278 {
279 bool badded = false;
280 Segment *seg;
281 node = m_seglist.GetFirst();
282 while (node)
283 {
284 seg = node->GetData();
285
286 if((!seg->bUsed) && (seg->py1 == tseg->py2) && (seg->px1 == tseg->px2)) // fits without reverse
287 {
288 seg->bUsed = true;
289 segjoin2.Append(seg);
290 badded = true;
291 break;
292 }
293 else if((!seg->bUsed) && (seg->py2 == tseg->py2) && (seg->px2 == tseg->px2)) // fits, needs reverse
294 {
295 seg->bUsed = true;
296 double a = seg->px2; seg->px2 = seg->px1; seg->px1 = a;
297 double b = seg->py2; seg->py2 = seg->py1; seg->py1 = b;
298 segjoin2.Append(seg);
299 badded = true;
300 break;
301 }
302
303 node = node->GetNext();
304 }
305 if(badded == true)
306 tseg = seg;
307 else
308 tseg = NULL;
309 }
310
311
312
313 // Build a chain extending from the "1" end of the target segment
314 // The joined list, side 1...
315 MySegList segjoin1;
316
317 // Add the same first segment to the list
318 node = m_seglist.GetFirst();
319 seg0 = node->GetData();
320 seg0->bUsed = true;
321 segjoin1.Append(seg0);
322
323 tseg = seg0;
324
325 while(tseg)
326 {
327 bool badded = false;
328 node = m_seglist.GetFirst();
329 while (node)
330 {
331 seg = node->GetData();
332
333 if((!seg->bUsed) && (seg->py2 == tseg->py1) && (seg->px2 == tseg->px1)) // fits without reverse
334 {
335 seg->bUsed = true;
336 segjoin1.Append(seg);
337 badded = true;
338 break;
339 }
340 else if((!seg->bUsed) && (seg->py1 == tseg->py1) && (seg->px1 == tseg->px1)) // fits, needs reverse
341 {
342 seg->bUsed = true;
343 double a = seg->px2; seg->px2 = seg->px1; seg->px1 = a;
344 double b = seg->py2; seg->py2 = seg->py1; seg->py1 = b;
345 segjoin1.Append(seg);
346 badded = true;
347 break;
348 }
349
350 node = node->GetNext();
351 }
352 if(badded == true)
353 tseg = seg;
354 else
355 tseg = NULL;
356 }
357
358
359 // Now have two lists...
360
361 // Start with "1" side list,
362 // starting from the end, and skipping the first segment
363
364 int n1 = segjoin1.GetCount();
365 for(int i=n1 - 1 ; i > 0 ; i--)
366 {
367 node = segjoin1.Item( i );
368 seg = node->GetData();
369 ret_list->Append(seg);
370 }
371
372 // Now add the "2"side list
373 int n2 = segjoin2.GetCount();
374 for(int i=0 ; i < n2 ; i++)
375 {
376 node = segjoin2.Item(i);
377 seg = node->GetData();
378 ret_list->Append(seg);
379 }
380
381 // And there it is
382
383 return ret_list;
384 }
385
386
387 //---------------------------------------------------------------
drawIsoLine(GRIBOverlayFactory * pof,wxDC * dc,PlugIn_ViewPort * vp,bool bHiDef)388 void IsoLine::drawIsoLine(GRIBOverlayFactory *pof, wxDC *dc, PlugIn_ViewPort *vp, bool bHiDef)
389 {
390 int nsegs = trace.size();
391 if(nsegs < 1)
392 return;
393
394 GetGlobalColor ( _T ( "UITX1" ), &isoLineColor );
395
396 #if wxUSE_GRAPHICS_CONTEXT
397 wxGraphicsContext *pgc = NULL;
398 #endif
399
400 if(dc) {
401 wxPen ppISO ( isoLineColor, 2 );
402
403 #if wxUSE_GRAPHICS_CONTEXT
404 wxMemoryDC *pmdc;
405 pmdc= wxDynamicCast(dc, wxMemoryDC);
406 pgc = wxGraphicsContext::Create(*pmdc);
407 pgc->SetPen(ppISO);
408 #endif
409 dc->SetPen(ppISO);
410 } else { /* opengl */
411 #ifdef ocpnUSE_GL
412 //#ifndef USE_ANDROID_GLES2
413 // if(m_pixelMM > 0.2){ // pixel size large enough to render well
414 // // Enable anti-aliased lines, at best quality
415 // glEnable( GL_LINE_SMOOTH );
416 // glEnable( GL_BLEND );
417 // glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
418 // glHint( GL_LINE_SMOOTH_HINT, GL_NICEST );
419 // glLineWidth( 2 );
420 // }
421 // else{
422 // glLineWidth( 0.4/m_pixelMM); // set a target line width by MM
423 // }
424 //#else
425 if(pof->m_oDC){
426 wxPen ppISO ( isoLineColor, 2 );
427 pof->m_oDC->SetPen(ppISO);
428 }
429 #endif
430 }
431
432 std::list<Segment *>::iterator it;
433
434 //---------------------------------------------------------
435 // Dessine les segments
436 //---------------------------------------------------------
437 for (it=trace.begin(); it!=trace.end(); it++)
438 {
439 Segment *seg = *it;
440
441 if(vp->m_projection_type == PI_PROJECTION_MERCATOR ||
442 vp->m_projection_type == PI_PROJECTION_EQUIRECTANGULAR) {
443 /* skip segments that go the wrong way around the world */
444 double sx1 = seg->px1, sx2 = seg->px2;
445 if(sx2 - sx1 > 180)
446 sx2 -= 360;
447 else if(sx1 - sx2 > 180)
448 sx1 -= 360;
449
450 if((sx1+180 < vp->clon && sx2+180 > vp->clon) ||
451 (sx1+180 > vp->clon && sx2+180 < vp->clon) ||
452 (sx1-180 < vp->clon && sx2-180 > vp->clon) ||
453 (sx1-180 > vp->clon && sx2-180 < vp->clon))
454 continue;
455 }
456
457 wxPoint ab;
458 GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
459 wxPoint cd;
460 GetCanvasPixLL(vp, &cd, seg->py2, seg->px2);
461
462 if(dc) {
463 #if wxUSE_GRAPHICS_CONTEXT
464 if(bHiDef && pgc)
465 pgc->StrokeLine(ab.x, ab.y, cd.x, cd.y);
466 else
467 #endif
468 dc->DrawLine(ab.x, ab.y, cd.x, cd.y);
469 } else { /* opengl */
470 #ifdef ocpnUSE_GL
471
472 if(pof->m_oDC){
473 pof->m_oDC->DrawLine(ab.x, ab.y, cd.x, cd.y);
474 }
475
476 #endif
477 }
478 }
479
480 #if wxUSE_GRAPHICS_CONTEXT
481 delete pgc;
482 #endif
483
484 // if(!dc) /* opengl */
485 // glEnd();
486 }
487
488 //---------------------------------------------------------------
489
drawIsoLineLabels(GRIBOverlayFactory * pof,wxDC * dc,PlugIn_ViewPort * vp,int density,int first,wxImage & imageLabel)490 void IsoLine::drawIsoLineLabels(GRIBOverlayFactory *pof, wxDC *dc,
491 PlugIn_ViewPort *vp, int density, int first,
492 wxImage &imageLabel)
493
494 {
495 std::list<Segment *>::iterator it;
496 int nb = first;
497 wxString label;
498
499 //---------------------------------------------------------
500 // Ecrit les labels
501 //---------------------------------------------------------
502 wxRect prev;
503 for (it=trace.begin(); it!=trace.end(); it++,nb++)
504 {
505 if (nb % density == 0)
506 {
507 Segment *seg = *it;
508
509 // if(vp->vpBBox.PointInBox((seg->px1 + seg->px2)/2., (seg->py1 + seg->py2)/2., 0.))
510 {
511 wxPoint ab;
512 GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
513 wxPoint cd;
514 GetCanvasPixLL(vp, &cd, seg->py1, seg->px1);
515
516 int w = imageLabel.GetWidth();
517 int h = imageLabel.GetHeight();
518
519 int label_offset = 6;
520 int xd = (ab.x + cd.x-(w+label_offset * 2))/2;
521 int yd = (ab.y + cd.y - h)/2;
522
523 int x = xd - label_offset;
524 wxRect r(x ,yd ,w ,h);
525 r.Inflate(w);
526 if (!prev.Intersects(r)) {
527 prev = r;
528
529 /* don't use alpha for isobars, for some reason draw bitmap ignores
530 the 4th argument (true or false has same result) */
531 wxImage img(w, h, imageLabel.GetData(), true);
532 dc->DrawBitmap(img, xd, yd, false);
533 }
534 }
535 }
536 }
537 }
538
drawIsoLineLabelsGL(GRIBOverlayFactory * pof,PlugIn_ViewPort * vp,int density,int first,wxString label,wxColour & color,TexFont & texfont)539 void IsoLine::drawIsoLineLabelsGL(GRIBOverlayFactory *pof,
540 PlugIn_ViewPort *vp, int density, int first,
541 wxString label, wxColour &color, TexFont &texfont)
542
543 {
544 std::list<Segment *>::iterator it;
545 int nb = first;
546
547 #ifdef ocpnUSE_GL
548 glEnable( GL_BLEND );
549 glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
550
551
552 //---------------------------------------------------------
553 // Ecrit les labels
554 //---------------------------------------------------------
555 wxRect prev;
556 for (it=trace.begin(); it!=trace.end(); it++,nb++)
557 {
558 if (nb % density == 0)
559 {
560 Segment *seg = *it;
561
562 // if(vp->vpBBox.PointInBox((seg->px1 + seg->px2)/2., (seg->py1 + seg->py2)/2., 0.))
563 {
564 wxPoint ab;
565 GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
566 wxPoint cd;
567 GetCanvasPixLL(vp, &cd, seg->py1, seg->px1);
568
569 int w, h;
570 texfont.GetTextExtent(label, &w, &h);
571
572 int label_offsetx = 6, label_offsety = 1;
573 int xd = (ab.x + cd.x-(w+label_offsetx * 2))/2;
574 int yd = (ab.y + cd.y - h)/2;
575 int x = xd - label_offsetx, y = yd - label_offsety;
576 w += 2*label_offsetx, h += 2*label_offsety;
577
578 wxRect r(x,y, w,h);
579 r.Inflate(w);
580 if (!prev.Intersects(r))
581 {
582 #if 1
583 prev = r;
584 if(pof->m_oDC){
585
586 //m_oDC->SetFont( *mfont );
587 pof->m_oDC->SetPen( *wxBLACK_PEN );
588 pof->m_oDC->SetBrush( color );
589 pof->m_oDC->DrawRectangle( x, y, w, h );
590 pof->m_oDC->DrawText( label, xd, yd );
591 }
592
593 #else
594 prev = r;
595 glColor4ub(color.Red(), color.Green(), color.Blue(), color.Alpha());
596
597 /* draw bounding rectangle */
598 glBegin(GL_QUADS);
599 glVertex2i(x, y);
600 glVertex2i(x+w, y);
601 glVertex2i(x+w, y+h);
602 glVertex2i(x, y+h);
603 glEnd();
604
605 glColor3ub(0, 0, 0);
606
607 glBegin(GL_LINE_LOOP);
608 glVertex2i(x, y);
609 glVertex2i(x+w, y);
610 glVertex2i(x+w, y+h);
611 glVertex2i(x, y+h);
612 glEnd();
613
614 glEnable(GL_TEXTURE_2D);
615 texfont.RenderString(label, xd, yd);
616 glDisable(GL_TEXTURE_2D);
617 #endif
618 }
619 }
620 }
621 }
622 glDisable( GL_BLEND );
623 #endif
624 }
625
626
627 //==================================================================================
628 // Segment
629 //==================================================================================
Segment(int I,int w,int J,char c1,char c2,char c3,char c4,const GribRecord * rec,double pressure)630 Segment::Segment(int I, int w, int J,
631 char c1, char c2, char c3, char c4,
632 const GribRecord *rec, double pressure)
633 {
634 traduitCode(I, w, J, c1, i,j);
635 traduitCode(I, w, J, c2, k,l);
636 traduitCode(I, w, J, c3, m,n);
637 traduitCode(I, w, J, c4, o,p);
638
639 intersectionAreteGrille(i,j, k,l, &px1,&py1, rec, pressure);
640 intersectionAreteGrille(m,n, o,p, &px2,&py2, rec, pressure);
641 }
642 //-----------------------------------------------------------------------
intersectionAreteGrille(int i,int j,int k,int l,double * x,double * y,const GribRecord * rec,double pressure)643 void Segment::intersectionAreteGrille(int i,int j, int k,int l, double *x, double *y,
644 const GribRecord *rec, double pressure)
645 {
646 double xa, xb, ya, yb, pa, pb, dec;
647 pa = rec->getValue(i,j);
648 pb = rec->getValue(k,l);
649
650 rec->getXY(i, j, &xa, &ya);
651 rec->getXY(k, l, &xb, &yb);
652
653 // Abscisse
654 if (pb != pa)
655 dec = (pressure-pa)/(pb-pa);
656 else
657 dec = 0.5;
658 if (fabs(dec)>1)
659 dec = 0.5;
660 double xd = xb -xa;
661 if(xd < -180)
662 xd += 360;
663 else if(xd > 180)
664 xd -= 360;
665 *x = xa +xd*dec;
666
667 // Ordonnée
668 if (pb != pa)
669 dec = (pressure-pa)/(pb-pa);
670 else
671 dec = 0.5;
672 if (fabs(dec)>1)
673 dec = 0.5;
674 *y = ya +(yb -ya)*dec;
675 }
676 //---------------------------------------------------------------
traduitCode(int I,int w,int J,char c1,int & i,int & j)677 void Segment::traduitCode(int I, int w, int J, char c1, int &i, int &j) {
678 int Im1 = I ? I-1 : w - 1;
679 switch (c1) {
680 case 'a': i=Im1; j=J-1; break;
681 case 'b': i=I ; j=J-1; break;
682 case 'c': i=Im1; j=J ; break;
683 case 'd': i=I ; j=J ; break;
684 default: i=I ; j=J ;
685 }
686 }
687
688 //-----------------------------------------------------------------------
689 // Génère la liste des segments.
690 // Les coordonnées sont les indices dans la grille du GribRecord
691 //---------------------------------------------------------
extractIsoLine(const GribRecord * rec)692 void IsoLine::extractIsoLine(const GribRecord *rec)
693 {
694 int i, j, W, H;
695 double a,b,c,d;
696 W = rec->getNi();
697 H = rec->getNj();
698
699 int We = W;
700 if(rec->getLonMax() + rec->getDi() - rec->getLonMin() == 360)
701 We++;
702
703 for (j=1; j<H; j++) // !!!! 1 to end
704 {
705 a = rec->getValue( 0, j-1 );
706 c = rec->getValue( 0, j );
707 for (i=1; i<We; i++, a = b, c = d)
708 {
709 // x = rec->getX(i);
710 // y = rec->getY(j);
711
712 int ni = i;
713 if (i == W)
714 ni = 0;
715 b = rec->getValue( ni, j-1 );
716 d = rec->getValue( ni, j );
717
718 if( a == GRIB_NOTDEF || b == GRIB_NOTDEF || c == GRIB_NOTDEF || d == GRIB_NOTDEF ) continue;
719
720 if ((a< value && b< value && c< value && d < value)
721 || (a>value && b>value && c>value && d > value))
722 continue;
723
724 // Détermine si 1 ou 2 segments traversent la case ab-cd
725 // a b
726 // c d
727 //--------------------------------
728 // 1 segment en diagonale
729 //--------------------------------
730 if ((a<=value && b<=value && c<=value && d>value)
731 || (a>value && b>value && c>value && d<=value))
732 trace.push_back(new Segment(ni,W,j, 'c','d', 'b','d', rec, value));
733 else if ((a<=value && c<=value && d<=value && b>value)
734 || (a>value && c>value && d>value && b<=value))
735 trace.push_back(new Segment(ni,W,j, 'a','b', 'b','d', rec, value));
736 else if ((c<=value && d<=value && b<=value && a>value)
737 || (c>value && d>value && b>value && a<=value))
738 trace.push_back(new Segment(ni,W,j, 'a','b', 'a','c', rec, value));
739 else if ((a<=value && b<=value && d<=value && c>value)
740 || (a>value && b>value && d>value && c<=value))
741 trace.push_back(new Segment(ni,W,j, 'a','c', 'c','d', rec,value));
742 //--------------------------------
743 // 1 segment H ou V
744 //--------------------------------
745 else if ((a<=value && b<=value && c>value && d>value)
746 || (a>value && b>value && c<=value && d<=value))
747 trace.push_back(new Segment(ni,W,j, 'a','c', 'b','d', rec,value));
748 else if ((a<=value && c<=value && b>value && d>value)
749 || (a>value && c>value && b<=value && d<=value))
750 trace.push_back(new Segment(ni,W,j, 'a','b', 'c','d', rec,value));
751 //--------------------------------
752 // 2 segments en diagonale
753 //--------------------------------
754 else if (a<=value && d<=value && c>value && b>value) {
755 trace.push_back(new Segment(ni,W,j, 'a','b', 'b','d', rec,value));
756 trace.push_back(new Segment(ni,W,j, 'a','c', 'c','d', rec,value));
757 }
758 else if (a>value && d>value && c<=value && b<=value) {
759 trace.push_back(new Segment(ni,W,j, 'a','b', 'a','c', rec,value));
760 trace.push_back(new Segment(ni,W,j, 'b','d', 'c','d', rec,value));
761 }
762
763 }
764 }
765 }
766
767
768 // ----------------------------------------------------------------------------
769 // splines code lifted from wxWidgets
770 // ----------------------------------------------------------------------------
771
772
773
774
775
776
777 // ----------------------------------- spline code ----------------------------------------
778
779 void ocpn_wx_quadratic_spline(double a1, double b1, double a2, double b2,
780 double a3, double b3, double a4, double b4);
781 void ocpn_wx_clear_stack();
782 int ocpn_wx_spline_pop(double *x1, double *y1, double *x2, double *y2, double *x3,
783 double *y3, double *x4, double *y4);
784 void ocpn_wx_spline_push(double x1, double y1, double x2, double y2, double x3, double y3,
785 double x4, double y4);
786 static bool ocpn_wx_spline_add_point(double x, double y);
787
788
789 #define half(z1, z2) ((z1+z2)/2.0)
790 #define THRESHOLD 5
791
792 /* iterative version */
793
ocpn_wx_quadratic_spline(double a1,double b1,double a2,double b2,double a3,double b3,double a4,double b4)794 void ocpn_wx_quadratic_spline(double a1, double b1, double a2, double b2, double a3, double b3, double a4,
795 double b4)
796 {
797 register double xmid, ymid;
798 double x1, y1, x2, y2, x3, y3, x4, y4;
799
800 ocpn_wx_clear_stack();
801 ocpn_wx_spline_push(a1, b1, a2, b2, a3, b3, a4, b4);
802
803 while (ocpn_wx_spline_pop(&x1, &y1, &x2, &y2, &x3, &y3, &x4, &y4)) {
804 xmid = (double)half(x2, x3);
805 ymid = (double)half(y2, y3);
806 if (fabs(x1 - xmid) < THRESHOLD && fabs(y1 - ymid) < THRESHOLD &&
807 fabs(xmid - x4) < THRESHOLD && fabs(ymid - y4) < THRESHOLD) {
808 ocpn_wx_spline_add_point( x1, y1 );
809 ocpn_wx_spline_add_point( xmid, ymid );
810 } else {
811 ocpn_wx_spline_push(xmid, ymid, (double)half(xmid, x3), (double)half(ymid, y3),
812 (double)half(x3, x4), (double)half(y3, y4), x4, y4);
813 ocpn_wx_spline_push(x1, y1, (double)half(x1, x2), (double)half(y1, y2),
814 (double)half(x2, xmid), (double)half(y2, ymid), xmid, ymid);
815 }
816 }
817 }
818
819 /* utilities used by spline drawing routines */
820
821 typedef struct ocpn_wx_spline_stack_struct {
822 double x1, y1, x2, y2, x3, y3, x4, y4;
823 } Stack;
824
825 #define SPLINE_STACK_DEPTH 20
826 static Stack ocpn_wx_spline_stack[SPLINE_STACK_DEPTH];
827 static Stack *ocpn_wx_stack_top;
828 static int ocpn_wx_stack_count;
829
ocpn_wx_clear_stack()830 void ocpn_wx_clear_stack()
831 {
832 ocpn_wx_stack_top = ocpn_wx_spline_stack;
833 ocpn_wx_stack_count = 0;
834 }
835
ocpn_wx_spline_push(double x1,double y1,double x2,double y2,double x3,double y3,double x4,double y4)836 void ocpn_wx_spline_push(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4)
837 {
838 ocpn_wx_stack_top->x1 = x1;
839 ocpn_wx_stack_top->y1 = y1;
840 ocpn_wx_stack_top->x2 = x2;
841 ocpn_wx_stack_top->y2 = y2;
842 ocpn_wx_stack_top->x3 = x3;
843 ocpn_wx_stack_top->y3 = y3;
844 ocpn_wx_stack_top->x4 = x4;
845 ocpn_wx_stack_top->y4 = y4;
846 ocpn_wx_stack_top++;
847 ocpn_wx_stack_count++;
848 }
849
ocpn_wx_spline_pop(double * x1,double * y1,double * x2,double * y2,double * x3,double * y3,double * x4,double * y4)850 int ocpn_wx_spline_pop(double *x1, double *y1, double *x2, double *y2,
851 double *x3, double *y3, double *x4, double *y4)
852 {
853 if (ocpn_wx_stack_count == 0)
854 return (0);
855 ocpn_wx_stack_top--;
856 ocpn_wx_stack_count--;
857 *x1 = ocpn_wx_stack_top->x1;
858 *y1 = ocpn_wx_stack_top->y1;
859 *x2 = ocpn_wx_stack_top->x2;
860 *y2 = ocpn_wx_stack_top->y2;
861 *x3 = ocpn_wx_stack_top->x3;
862 *y3 = ocpn_wx_stack_top->y3;
863 *x4 = ocpn_wx_stack_top->x4;
864 *y4 = ocpn_wx_stack_top->y4;
865 return (1);
866 }
867
ocpn_wx_spline_add_point(double x,double y)868 static bool ocpn_wx_spline_add_point(double x, double y)
869 {
870 wxPoint *point = new wxPoint ;
871 point->x = (int) x;
872 point->y = (int) y;
873 ocpn_wx_spline_point_list.Append((wxObject*)point);
874 return true;
875 }
876
877
GenSpline(wxList * points)878 void GenSpline( wxList *points )
879 {
880
881 wxPoint *p;
882 double cx1, cy1, cx2, cy2, cx3, cy3, cx4, cy4;
883 double x1, y1, x2, y2;
884
885 wxList::compatibility_iterator node = points->GetFirst();
886 if (!node)
887 // empty list
888 return;
889
890 p = (wxPoint *)node->GetData();
891
892 x1 = p->x;
893 y1 = p->y;
894
895 node = node->GetNext();
896 p = (wxPoint *)node->GetData();
897
898 x2 = p->x;
899 y2 = p->y;
900 cx1 = (double)((x1 + x2) / 2);
901 cy1 = (double)((y1 + y2) / 2);
902 cx2 = (double)((cx1 + x2) / 2);
903 cy2 = (double)((cy1 + y2) / 2);
904
905 ocpn_wx_spline_add_point(x1, y1);
906
907 while ((node = node->GetNext())
908 #if !wxUSE_STL
909 != NULL
910 #endif // !wxUSE_STL
911 )
912 {
913 p = (wxPoint *)node->GetData();
914 x1 = x2;
915 y1 = y2;
916 x2 = p->x;
917 y2 = p->y;
918 cx4 = (double)(x1 + x2) / 2;
919 cy4 = (double)(y1 + y2) / 2;
920 cx3 = (double)(x1 + cx4) / 2;
921 cy3 = (double)(y1 + cy4) / 2;
922
923 ocpn_wx_quadratic_spline(cx1, cy1, cx2, cy2, cx3, cy3, cx4, cy4);
924
925 cx1 = cx4;
926 cy1 = cy4;
927 cx2 = (double)(cx1 + x2) / 2;
928 cy2 = (double)(cy1 + y2) / 2;
929 }
930
931 ocpn_wx_spline_add_point( cx1, cy1 );
932 ocpn_wx_spline_add_point( x2, y2 );
933
934 }
935
936 #if 0
937 static void GenerateSpline(int n, wxPoint points[])
938 {
939 wxList list;
940 for (int i =0; i < n; i++)
941 {
942 list.Append((wxObject*)&points[i]);
943 }
944
945 GenSpline(&list);
946 }
947
948 static void ClearSplineList()
949 {
950 wxList::compatibility_iterator node = ocpn_wx_spline_point_list.GetFirst();
951 while (node)
952 {
953 wxPoint *point = (wxPoint *)node->GetData();
954 delete point;
955 ocpn_wx_spline_point_list.Erase(node);
956 node = ocpn_wx_spline_point_list.GetFirst();
957 }
958 }
959 #endif
960