@@ -191,4 +191,116 @@ void InverseViewshed::addEventsFromCell( int &row, int &column, const double &pi
191191 mCellEvents .push_back ( mEventCenter );
192192 }
193193 }
194+ }
195+
196+ void InverseViewshed::calculateVisibilityRaster ()
197+ {
198+
199+ mVisibilityRaster = std::make_shared<ProjectedSquareCellRaster>( *mInputDsm .get (), GDALDataType::GDT_Float32 );
200+
201+ std::chrono::_V2::steady_clock::time_point startTime = std::chrono::_V2::steady_clock::now ();
202+
203+ std::map<double , LoSNode> los;
204+
205+ for ( LoSNode ln : mLosNodes )
206+ {
207+ if ( ln.mInverseLoSBehindTarget )
208+ {
209+ continue ;
210+ }
211+
212+ los.insert ( std::pair<double , LoSNode>( ln.centreDistance (), ln ) );
213+ }
214+
215+ std::size_t i = 0 ;
216+ for ( CellEvent e : mCellEvents )
217+ {
218+ // progressCallback( mCellEvents.size(), i );
219+
220+ switch ( e.mEventType )
221+ {
222+ case CellEventPositionType::ENTER :
223+ {
224+ if ( mPoint ->mRow == e.mRow && mPoint ->mCol == e.mCol )
225+ {
226+ break ;
227+ }
228+ if ( e.mBehindTargetForInverseLoS )
229+ {
230+ break ;
231+ }
232+
233+ mLoSNodeTemp = LoSNode ( mPoint ->mRow , mPoint ->mCol , &e, mCellSize );
234+ los.insert ( std::pair<double , LoSNode>( mLoSNodeTemp .centreDistance (), mLoSNodeTemp ) );
235+ break ;
236+ }
237+ case CellEventPositionType::EXIT :
238+ {
239+ if ( mPoint ->mRow == e.mRow && mPoint ->mCol == e.mCol )
240+ {
241+ break ;
242+ }
243+ if ( e.mBehindTargetForInverseLoS )
244+ {
245+ break ;
246+ }
247+
248+ mLoSNodeTemp = LoSNode ( mPoint ->mRow , mPoint ->mCol , &e, mCellSize );
249+ los.erase ( mLoSNodeTemp .centreDistance () );
250+ break ;
251+ }
252+ case CellEventPositionType::CENTER :
253+ {
254+ double maxGradient = -99 ;
255+
256+ mLoSNodeTemp = LoSNode ( mPoint ->mRow , mPoint ->mCol , &e, mCellSize );
257+ double curCorr, gradient, elevation, horizontalAngle;
258+
259+ double visibilityElevation = mLoSNodeTemp .centreElevation () + mObserverOffset ;
260+
261+ curCorr = Visibility::curvatureCorrections ( mLoSNodeTemp .centreDistance (), mRefractionCoefficient ,
262+ mEarthDiameter );
263+ elevation = visibilityElevation + curCorr;
264+
265+ gradient = Visibility::gradient ( mPoint ->totalElevation () - elevation, mLoSNodeTemp .centreDistance () );
266+
267+ double visibilityGradient = gradient;
268+ double visibilityDistance = mLoSNodeTemp .centreDistance ();
269+
270+ for ( auto it = los.begin (); it != los.lower_bound ( mLoSNodeTemp .centreDistance () ); ++it )
271+ {
272+ horizontalAngle = mLoSNodeTemp .centreAngle (); // + M_PI;
273+
274+ curCorr = Visibility::curvatureCorrections ( visibilityDistance -
275+ it->second .distanceAtAngle ( horizontalAngle ),
276+ mRefractionCoefficient , mEarthDiameter );
277+
278+ elevation = it->second .elevationAtAngle ( horizontalAngle ) + curCorr;
279+
280+ gradient =
281+ Visibility::gradient ( elevation - visibilityElevation,
282+ visibilityDistance - it->second .distanceAtAngle ( horizontalAngle ) );
283+
284+ if ( gradient > maxGradient )
285+ {
286+ maxGradient = gradient;
287+ }
288+ }
289+
290+ double visible = 0 ;
291+ if ( maxGradient < visibilityGradient )
292+ {
293+ visible = 1 ;
294+ };
295+ mVisibilityRaster ->writeValue ( e.mRow , e.mCol , visible );
296+
297+ break ;
298+ }
299+ }
300+
301+ i++;
302+ }
303+
304+ mTimeParse =
305+ std::chrono::duration_cast<std::chrono::nanoseconds>( std::chrono::_V2::steady_clock::now () - startTime );
194306}
0 commit comments