1 /* -*-c++-*- */
2 /* osgEarth - Geospatial SDK for OpenSceneGraph
3  * Copyright 2019 Pelican Mapping
4  * http://osgearth.org
5  *
6  * osgEarth is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with this program.  If not, see <http://www.gnu.org/licenses/>
18  */
19 
20 #include <osgEarthFeatures/FeatureTileSource>
21 #include <osgEarthFeatures/ResampleFilter>
22 #include <osgEarthFeatures/TransformFilter>
23 #include <osgEarthFeatures/BufferFilter>
24 #include <osgEarthFeatures/FilterContext>
25 #include <osgEarthSymbology/Style>
26 //TODO: replace this with GeometryRasterizer
27 #include <osgEarthSymbology/AGG.h>
28 #include <osgEarth/Registry>
29 #include <osgEarth/FileUtils>
30 #include <osgEarth/ImageUtils>
31 
32 #include <osg/Notify>
33 #include <osgDB/FileNameUtils>
34 #include <osgDB/FileUtils>
35 #include <osgDB/Registry>
36 #include <osgDB/ReadFile>
37 #include <osgDB/WriteFile>
38 
39 #include "AGGLiteOptions"
40 
41 #include <sstream>
42 #include <OpenThreads/Mutex>
43 #include <OpenThreads/ScopedLock>
44 
45 #define LC "[AGGLite] "
46 
47 using namespace osgEarth;
48 using namespace osgEarth::Features;
49 using namespace osgEarth::Symbology;
50 using namespace osgEarth::Drivers;
51 using namespace OpenThreads;
52 
53 namespace
54 {
55     struct float32
56     {
float32__anond9f1513f0111::float3257         float32() : value(NO_DATA_VALUE) { }
float32__anond9f1513f0111::float3258         float32(float v) : value(v) { }
59 
60         float value;
61     };
62 
63     struct span_coverage32
64     {
render__anond9f1513f0111::span_coverage3265         static void render(unsigned char* ptr,
66                            int x,
67                            unsigned count,
68                            const unsigned char* covers,
69                            const float32& c)
70         {
71             unsigned char* p = ptr + (x << 2);
72             float* f = (float*)p;
73             do
74             {
75                 unsigned char cover = *covers++;
76                 int hasData = cover > 127;
77                 *f++ = hasData ? c.value : NO_DATA_VALUE;
78             }
79             while(--count);
80         }
81 
hline__anond9f1513f0111::span_coverage3282         static void hline(unsigned char* ptr,
83                           int x,
84                           unsigned count,
85                           const float32& c)
86         {
87             unsigned char* p = ptr + (x << 2);
88             float* f = (float*)p;
89             do {
90                 *f++ = c.value;
91             }
92             while(--count);
93         }
94 
get__anond9f1513f0111::span_coverage3295         static float32 get(unsigned char* ptr, int x)
96         {
97             unsigned char* p = ptr + (x << 2);
98             float* f = (float*)p;
99             return float32(*f);
100         }
101     };
102 }
103 
104 /********************************************************************/
105 
106 class AGGLiteRasterizerTileSource : public FeatureTileSource
107 {
108 public:
109     struct RenderFrame {
110         double xmin, ymin;
111         double xf, yf;
112     };
113 
114 public:
AGGLiteRasterizerTileSource(const TileSourceOptions & options)115     AGGLiteRasterizerTileSource( const TileSourceOptions& options ) : FeatureTileSource( options ),
116         _options( options )
117     {
118         //nop
119     }
120 
121     //override
allocateImage()122     osg::Image* allocateImage()
123     {
124         osg::Image* image = 0L;
125         if ( _options.coverage() == true )
126         {
127             image = new osg::Image();
128             image->allocateImage(getPixelsPerTile(), getPixelsPerTile(), 1, GL_LUMINANCE, GL_FLOAT);
129             image->setInternalTextureFormat(GL_LUMINANCE32F_ARB);
130             ImageUtils::markAsUnNormalized(image, true);
131         }
132         return image;
133     }
134 
135     //override
preProcess(osg::Image * image,osg::Referenced * buildData)136     bool preProcess(osg::Image* image, osg::Referenced* buildData)
137     {
138         agg::rendering_buffer rbuf( image->data(), image->s(), image->t(), image->s()*4 );
139 
140         // clear the buffer.
141         if ( _options.coverage() == true )
142         {
143             agg::renderer<span_coverage32, float32> ren(rbuf);
144             ren.clear( float32(NO_DATA_VALUE) );
145         }
146         else
147         {
148             agg::renderer<agg::span_abgr32, agg::rgba8> ren(rbuf);
149             ren.clear(agg::rgba8(0,0,0,0));
150         }
151         return true;
152     }
153 
154     //override
renderFeaturesForStyle(Session * session,const Style & style,const FeatureList & features,osg::Referenced * buildData,const GeoExtent & imageExtent,osg::Image * image)155     bool renderFeaturesForStyle(
156         Session*           session,
157         const Style&       style,
158         const FeatureList& features,
159         osg::Referenced*   buildData,
160         const GeoExtent&   imageExtent,
161         osg::Image*        image )
162     {
163         OE_DEBUG << LC << "Rendering " << features.size() << " features for " << imageExtent.toString() << "\n";
164 
165         // A processing context to use with the filters:
166         FilterContext context( session );
167         context.setProfile( getFeatureSource()->getFeatureProfile() );
168 
169         const LineSymbol*    masterLine = style.getSymbol<LineSymbol>();
170         const PolygonSymbol* masterPoly = style.getSymbol<PolygonSymbol>();
171         const CoverageSymbol* masterCov = style.getSymbol<CoverageSymbol>();
172 
173         // sort into bins, making a copy for lines that require buffering.
174         FeatureList polygons;
175         FeatureList lines;
176 
177         for(FeatureList::const_iterator f = features.begin(); f != features.end(); ++f)
178         {
179             if ( f->get()->getGeometry() )
180             {
181                 bool hasPoly = false;
182                 bool hasLine = false;
183 
184                 if ( masterPoly || f->get()->style()->has<PolygonSymbol>() )
185                 {
186                     polygons.push_back( f->get() );
187                     hasPoly = true;
188                 }
189 
190                 if ( masterLine || f->get()->style()->has<LineSymbol>() )
191                 {
192                     // Use the GeometryIterator to get all the geometries so we can clone them as rings
193                     GeometryIterator gi(f->get()->getGeometry());
194                     while (gi.hasMore())
195                     {
196                         Geometry* geom = gi.next();
197                         // Create a new feature for each geometry
198                         Feature* newFeature = new Feature(*f->get());
199                         newFeature->setGeometry(geom);
200                         if (!newFeature->getGeometry()->isLinear())
201                         {
202                             newFeature->setGeometry(newFeature->getGeometry()->cloneAs(Geometry::TYPE_RING));
203                         }
204                         lines.push_back( newFeature );
205                         hasLine = true;
206                     }
207                 }
208 
209                 // if there are no geometry symbols but there is a coverage symbol, default to polygons.
210                 if ( !hasLine && !hasPoly )
211                 {
212                     if ( masterCov || f->get()->style()->has<CoverageSymbol>() )
213                     {
214                         polygons.push_back( f->get() );
215                     }
216                 }
217             }
218         }
219 
220         // initialize:
221         RenderFrame frame;
222         frame.xmin = imageExtent.xMin();
223         frame.ymin = imageExtent.yMin();
224         frame.xf   = (double)image->s() / imageExtent.width();
225         frame.yf   = (double)image->t() / imageExtent.height();
226 
227         if ( lines.size() > 0 )
228         {
229             // We are buffering in the features native extent, so we need to use the
230             // transformed extent to get the proper "resolution" for the image
231             const SpatialReference* featureSRS = context.profile()->getSRS();
232             GeoExtent transformedExtent = imageExtent.transform(featureSRS);
233 
234             double trans_xf = (double)image->s() / transformedExtent.width();
235             double trans_yf = (double)image->t() / transformedExtent.height();
236 
237             // resolution of the image (pixel extents):
238             double xres = 1.0/trans_xf;
239             double yres = 1.0/trans_yf;
240 
241             // downsample the line data so that it is no higher resolution than to image to which
242             // we intend to rasterize it. If you don't do this, you run the risk of the buffer
243             // operation taking forever on very high-res input data.
244             if ( _options.optimizeLineSampling() == true )
245             {
246                 ResampleFilter resample;
247                 resample.minLength() = osg::minimum( xres, yres );
248                 context = resample.push( lines, context );
249             }
250 
251             // now run the buffer operation on all lines:
252             BufferFilter buffer;
253             double lineWidth = 1.0;
254             if ( masterLine )
255             {
256                 buffer.capStyle() = masterLine->stroke()->lineCap().value();
257 
258                 if ( masterLine->stroke()->width().isSet() )
259                 {
260                     lineWidth = masterLine->stroke()->width().value();
261 
262                     GeoExtent imageExtentInFeatureSRS = imageExtent.transform(featureSRS);
263                     double pixelWidth = imageExtentInFeatureSRS.width() / (double)image->s();
264 
265                     // if the width units are specified, process them:
266                     if (masterLine->stroke()->widthUnits().isSet() &&
267                         masterLine->stroke()->widthUnits().get() != Units::PIXELS)
268                     {
269                         const Units& featureUnits = featureSRS->getUnits();
270                         const Units& strokeUnits  = masterLine->stroke()->widthUnits().value();
271 
272                         // if the units are different than those of the feature data, we need to
273                         // do a units conversion.
274                         if ( featureUnits != strokeUnits )
275                         {
276                             if ( Units::canConvert(strokeUnits, featureUnits) )
277                             {
278                                 // linear to linear, no problem
279                                 lineWidth = strokeUnits.convertTo( featureUnits, lineWidth );
280                             }
281                             else if ( strokeUnits.isLinear() && featureUnits.isAngular() )
282                             {
283                                 // linear to angular? approximate degrees per meter at the
284                                 // latitude of the tile's centroid.
285                                 double lineWidthM = masterLine->stroke()->widthUnits()->convertTo(Units::METERS, lineWidth);
286                                 double mPerDegAtEquatorInv = 360.0/(featureSRS->getEllipsoid()->getRadiusEquator() * 2.0 * osg::PI);
287                                 double lon, lat;
288                                 imageExtent.getCentroid(lon, lat);
289                                 lineWidth = lineWidthM * mPerDegAtEquatorInv * cos(osg::DegreesToRadians(lat));
290                             }
291                         }
292 
293                         // enfore a minimum width of one pixel.
294                         float minPixels = masterLine->stroke()->minPixels().getOrUse( 1.0f );
295                         lineWidth = osg::clampAbove(lineWidth, pixelWidth*minPixels);
296                     }
297 
298                     else // pixels
299                     {
300                         lineWidth *= pixelWidth;
301                     }
302                 }
303             }
304 
305             buffer.distance() = lineWidth * 0.5;   // since the distance is for one side
306             buffer.push( lines, context );
307         }
308 
309         // Transform the features into the map's SRS:
310         TransformFilter xform( imageExtent.getSRS() );
311         xform.setLocalizeCoordinates( false );
312         FilterContext polysContext = xform.push( polygons, context );
313         FilterContext linesContext = xform.push( lines, context );
314 
315         // set up the AGG renderer:
316         agg::rendering_buffer rbuf( image->data(), image->s(), image->t(), image->s()*4 );
317 
318         // Create the renderer and the rasterizer
319         agg::rasterizer ras;
320 
321         // Setup the rasterizer
322         if ( _options.coverage() == true )
323             ras.gamma(1.0);
324         else
325             ras.gamma(_options.gamma().get());
326 
327         ras.filling_rule(agg::fill_even_odd);
328 
329         // construct an extent for cropping the geometry to our tile.
330         // extend just outside the actual extents so we don't get edge artifacts:
331         GeoExtent cropExtent = GeoExtent(imageExtent);
332         cropExtent.scale(1.1, 1.1);
333         double cropXMin, cropYMin, cropXMax, cropYMax;
334         cropExtent.getBounds(cropXMin, cropYMin, cropXMax, cropYMax);
335 
336         // GEOS crop won't abide by weird extents, so if we're in geographic space
337         // we must clamp the scaled extent back to a legal range.
338         if (cropExtent.crossesAntimeridian())
339         {
340             osg::Vec3d centroid = imageExtent.getCentroid();
341             if (centroid.x() < 0.0) // tile is east of antimeridian
342             {
343                 cropXMin = -180.0;
344                 cropXMax = cropExtent.east();
345             }
346             else
347             {
348                 cropXMin = cropExtent.west();
349                 cropXMax = 180.0;
350             }
351         }
352 
353         osg::ref_ptr<Symbology::Polygon> cropPoly = new Symbology::Polygon( 4 );
354         cropPoly->push_back( osg::Vec3d(cropXMin, cropYMin, 0) );
355         cropPoly->push_back( osg::Vec3d(cropXMax, cropYMin, 0) );
356         cropPoly->push_back( osg::Vec3d(cropXMax, cropYMax, 0) );
357         cropPoly->push_back( osg::Vec3d(cropXMin, cropYMax, 0) );
358 
359         // If there's a coverage symbol, make a copy of the expressions so we can evaluate them
360         optional<NumericExpression> covValue;
361         const CoverageSymbol* covsym = style.get<CoverageSymbol>();
362         if (covsym && covsym->valueExpression().isSet())
363             covValue = covsym->valueExpression().get();
364 
365         // render the polygons
366         for(FeatureList::iterator i = polygons.begin(); i != polygons.end(); i++)
367         {
368             Feature*  feature  = i->get();
369             Geometry* geometry = feature->getGeometry();
370 
371             osg::ref_ptr<Geometry> croppedGeometry;
372             if ( geometry->crop( cropPoly.get(), croppedGeometry ) )
373             {
374                 const PolygonSymbol* poly =
375                     feature->style().isSet() && feature->style()->has<PolygonSymbol>() ? feature->style()->get<PolygonSymbol>() :
376                     masterPoly;
377 
378                 if ( _options.coverage() == true && covValue.isSet() )
379                 {
380                     float value = (float)feature->eval(covValue.mutable_value(), &context);
381                     rasterizeCoverage(croppedGeometry.get(), value, frame, ras, rbuf);
382                 }
383                 else
384                 {
385                     osg::Vec4f color = poly->fill()->color();
386                     rasterize(croppedGeometry.get(), color, frame, ras, rbuf);
387                 }
388 
389             }
390         }
391 
392         // render the lines
393         for(FeatureList::iterator i = lines.begin(); i != lines.end(); i++)
394         {
395             Feature*  feature  = i->get();
396             Geometry* geometry = feature->getGeometry();
397 
398             osg::ref_ptr<Geometry> croppedGeometry;
399             if ( geometry->crop( cropPoly.get(), croppedGeometry ) )
400             {
401                 const LineSymbol* line =
402                     feature->style().isSet() && feature->style()->has<LineSymbol>() ? feature->style()->get<LineSymbol>() :
403                     masterLine;
404 
405                 if ( _options.coverage() == true && covValue.isSet() )
406                 {
407                     float value = (float)feature->eval(covValue.mutable_value(), &context);
408                     rasterizeCoverage(croppedGeometry.get(), value, frame, ras, rbuf);
409                 }
410                 else
411                 {   osg::Vec4f color = line ? static_cast<osg::Vec4>(line->stroke()->color()) : osg::Vec4(1,1,1,1);
412                     rasterize(croppedGeometry.get(), color, frame, ras, rbuf);
413                 }
414             }
415         }
416 
417         return true;
418     }
419 
420     //override
postProcess(osg::Image * image,osg::Referenced * data)421     bool postProcess( osg::Image* image, osg::Referenced* data )
422     {
423         if ( _options.coverage() == false )
424         {
425             //convert from ABGR to RGBA
426             unsigned char* pixel = image->data();
427             for(int i=0; i<image->s()*image->t()*4; i+=4, pixel+=4)
428             {
429                 std::swap( pixel[0], pixel[3] );
430                 std::swap( pixel[1], pixel[2] );
431             }
432         }
433 
434         return true;
435     }
436 
437     // rasterizes a geometry to color
rasterize(const Geometry * geometry,const osg::Vec4 & color,RenderFrame & frame,agg::rasterizer & ras,agg::rendering_buffer & buffer)438     void rasterize(const Geometry* geometry, const osg::Vec4& color, RenderFrame& frame,
439                    agg::rasterizer& ras, agg::rendering_buffer& buffer)
440     {
441         unsigned a = (unsigned)(127.0f+(color.a()*255.0f)/2.0f); // scale alpha up
442         agg::rgba8 fgColor = agg::rgba8( (unsigned)(color.r()*255.0f), (unsigned)(color.g()*255.0f), (unsigned)(color.b()*255.0f), a );
443 
444         ConstGeometryIterator gi( geometry );
445         while( gi.hasMore() )
446         {
447             const Geometry* g = gi.next();
448 
449             for( Geometry::const_iterator p = g->begin(); p != g->end(); p++ )
450             {
451                 const osg::Vec3d& p0 = *p;
452                 double x0 = frame.xf*(p0.x()-frame.xmin);
453                 double y0 = frame.yf*(p0.y()-frame.ymin);
454 
455                 if ( p == g->begin() )
456                     ras.move_to_d( x0, y0 );
457                 else
458                     ras.line_to_d( x0, y0 );
459             }
460         }
461         agg::renderer<agg::span_abgr32, agg::rgba8> ren(buffer);
462         ras.render(ren, fgColor);
463 
464         ras.reset();
465     }
466 
467 
rasterizeCoverage(const Geometry * geometry,float value,RenderFrame & frame,agg::rasterizer & ras,agg::rendering_buffer & buffer)468     void rasterizeCoverage(const Geometry* geometry, float value, RenderFrame& frame,
469                            agg::rasterizer& ras, agg::rendering_buffer& buffer)
470     {
471         ConstGeometryIterator gi( geometry );
472         while( gi.hasMore() )
473         {
474             const Geometry* g = gi.next();
475 
476             for( Geometry::const_iterator p = g->begin(); p != g->end(); p++ )
477             {
478                 const osg::Vec3d& p0 = *p;
479                 double x0 = frame.xf*(p0.x()-frame.xmin);
480                 double y0 = frame.yf*(p0.y()-frame.ymin);
481 
482                 if ( p == g->begin() )
483                     ras.move_to_d( x0, y0 );
484                 else
485                     ras.line_to_d( x0, y0 );
486             }
487         }
488 
489         agg::renderer<span_coverage32, float32> ren(buffer);
490         ras.render(ren, value);
491         ras.reset();
492     }
493 
494 
getExtension() const495     virtual std::string getExtension()  const
496     {
497         return "png";
498     }
499 
500 private:
501     const AGGLiteOptions _options;
502     std::string _configPath;
503 };
504 
505 
506 /**
507  * Plugin entry point for the AGGLite feature rasterizer
508  */
509 class AGGLiteRasterizerTileSourceDriver : public TileSourceDriver
510 {
511     public:
AGGLiteRasterizerTileSourceDriver()512         AGGLiteRasterizerTileSourceDriver() {}
513 
className() const514         virtual const char* className() const
515         {
516             return "AGG-Lite feature rasterizer";
517         }
518 
acceptsExtension(const std::string & extension) const519         virtual bool acceptsExtension(const std::string& extension) const
520         {
521             return
522                 osgDB::equalCaseInsensitive( extension, "osgearth_agglite" ) ||
523                 osgDB::equalCaseInsensitive( extension, "osgearth_rasterize" );
524         }
525 
readObject(const std::string & file_name,const Options * options) const526         virtual ReadResult readObject(const std::string& file_name, const Options* options) const
527         {
528             std::string ext = osgDB::getFileExtension( file_name );
529             if ( !acceptsExtension( ext ) )
530             {
531                 return ReadResult::FILE_NOT_HANDLED;
532             }
533 
534             return new AGGLiteRasterizerTileSource( getTileSourceOptions(options) );
535         }
536 };
537 
538 REGISTER_OSGPLUGIN(osgearth_agglite, AGGLiteRasterizerTileSourceDriver)
539