ProductFQD.h

 
#ifndef POLARMAX_H_
#define POLARMAX_H_
 
#include "Geometry.h"  // radar geometry
 
#include "drain/image/Image.h"
#include "drain/image/ImageView.h"
#include "drain/image/ImageOp.h"
 
//#include "geometry.h"
 
namespace drain
{
 
namespace radar
{
 
 
    
template <class T = unsigned char,class T2 = unsigned char>
class ProductFQD : public image::ImageOp<T,T2>
{
public:
    ProductFQD(){
        this->setInfo("Computes maximum altitude of echo.","","");
    };
    //virtual ~ProductFQD(){};
    image::Image<T> &filter(const image::Image<T> &src,image::Image<T2> &dst) const;
};
 
template <class T,class T2>
image::Image<T> &ProductFQD<T,T2>::filter(const image::Image<T> &src,image::Image<T2> &dst) const
{
   const unsigned int srcWidth  = src.getWidth();
   const unsigned int srcHeight = src.getHeight();
   
   const unsigned int width =  this->parameters.get("width",src.getWidth());
   const unsigned int height = srcHeight; 
   dst.setGeometry(width,height,1,1);
   
   radar::Geometry radarGeometry;
   float binDepth = src.properties.get("BINDEPTH",500.0);
   
   Data elevations = src.properties.get("ELEVATIONS","0.123");
   elevations.trim();
   elevations.splitTo(radarGeometry.elevationAngles," ");
   
   const unsigned int sweeps = radarGeometry.elevationAngles.size();
   
   Data bins = src.properties.get("BINS");
   bins.trim();
   bins.splitTo(radarGeometry.bins," ");
   if (radarGeometry.bins.size() != sweeps){
        radarGeometry.bins.resize(sweeps,srcWidth);
   } 
   
   for (unsigned int i = 0; i < sweeps; ++i) {
        radarGeometry.elevationAngles[i] *= M_PI/180.0;
   }
    
   
    unsigned int bin;
 
    unsigned int groundBin;
 
    float eta;
 
    float beta;
 
    
    if (false){
        std::cout << "Starting Max\n";
        std::cout << "\nParams:\n"; 
        std::cout << this->parameters << '\n';
    
        std::cout << "\nGeometry:\n";   
        std::cout << src.getGeometry() << '\n';
        std::cout << "\nProperties:\n"; 
        std::cout << src.properties << '\n';
        std::cout << "\nElevs: " << elevations << "\n";
    
        for (unsigned int i=0; i< radarGeometry.elevationAngles.size(); i++){
            std::cout << "'" <<  radarGeometry.elevationAngles[i] << "', "; 
        }
        std::cout << "\n";
    }
    
    double groundDistance = 0.0;
    
    
    // Main loop
    for (unsigned int k=0; k<sweeps; k++){
 
      const drain::image::Image<T> & srcSweep = src.getChannel(k);
      const unsigned int srcWidth  = srcSweep.getWidth();
    
      eta  = radarGeometry.elevationAngles[k];
      
      for (unsigned int i=0; i<width; i++){
 
        // Scale in width
        groundBin = (i * srcWidth) / width;
      
        groundDistance = groundBin*binDepth + binDepth/2.0;
        beta = Geometry::betaFromGround(groundDistance);
        
        bin =  static_cast<unsigned int>(Geometry::beamFromEtaBeta(eta,beta) / binDepth);
        
        if (bin < radarGeometry.bins[k]) {
        
            for (unsigned int j=0; j<height; j++){
                T2  s = static_cast<T2>(srcSweep.at(bin,j));
                T2 &d = dst.at(bin,j);
                if (s > d) 
                  d = s;    
            }
        
        }
       
      }
      //std::cerr << "Altitude at end ("<< groundDistance<<"m): " << altitude/scale << std::endl;
    }
   
   
   return dst;  
};
 
 
} // ::image
 
} // ::drain
 
#endif /*POLARMAX_H_*/