AnoRack.h

 
#ifndef ANORACK_H_
#define ANORACK_H_ "AnoRack v0.1"
 
#include "radar/Constants.h"  // radar
#include "radar/Andre.h"
 
using namespace drain::image;
 
namespace anorack
{
 
 
template <class T = unsigned char,class T2 = unsigned char>
class AnoRackShip : public rack::AndreDetector<T,T2>
{
public:
 
 
    AnoRackShip(const std::string & p="64,4,5")
    : rack::AndreDetector<T,T2>("ship","Detects ships","threshold,size,sidelobe",p){
 
    };
 
 
    void filter(const Image<T> &src, Image<T2> &dst) const {
 
        makeCompatible(src,dst);
 
        const int threshold = this->getParameter("threshold",64);
        const int size = this->getParameter("size",4);
        const int sidelobe = this->getParameter("sidelobe",5);
 
        /*
            HighBoostOp<T> high;
            high.setParameter("width",3);
            high.setParameter("height",3);
            //high.filter(src,this->tmp);
            high.filter(src,dst);
         */
 
        SegmentAreaOp<T,T> area;
        area.setParameter("mapping","f");  // bilinear inverse
        area.setParameter("min",threshold);
        area.setParameter("targetSize",size);
        area.setParameter("halfWidth",size/2);
        area.filter(src,dst);
        MultiplicationOp<T,T>().filter(src,dst,dst);
 
        if (sidelobe > 0){
            DistanceTransformExponentialOp<T,T> dist;
            //DistanceTransformLinearOp<T,T> dist;
            dist.setParameter("horz",1);
            dist.setParameter("vert",sidelobe);
            dist.setParameter("diag",1);
            dist.filter(dst,dst);
        }
 
    };
};
 
 
template <class T = unsigned char,class T2 = unsigned char>
class AnoRackRadomeAttenuation : public rack::AndreDetector<T,T2>
{
public:
 
 
    AnoRackRadomeAttenuation(const std::string & p="45,10,20,2")
    : rack::AndreDetector<T,T2>("aRadome","Detects radome attenuation.","dBZmax,binStart,binEnd,power",p) {
    };
 
 
    void filter(const Image<T> &src,Image<T2> &dst) const {
 
        makeCompatible(src,dst);
 
        const unsigned int width = src.getWidth();
        const unsigned int height = src.getHeight();
 
        const float dBZmax = this->getParameter("dBZmax",45);
        const unsigned int binStart = std::min(this->getParameter("binStart",10u),width-1);
        const unsigned int binEnd   = std::min(this->getParameter("binEnd",20u),width-1);
        const double power  = this->getParameter("power",2.0);
 
        // TODO generalize FMI scale
        const double byteMax = dBZmax*2.0 + 64.0;
 
        std::cerr << binStart << ':' << binEnd << '\n';
 
        double sum = 0.0;
        double x;
        for (unsigned int i = binStart; i < binEnd; ++i) {
            for (unsigned int j = 0; j < height; ++j) {
                x = src.at(i,j);
                x = std::min(x,byteMax);
                //dst.at(i,j) = x;
                sum += ::pow(x/byteMax,power);
            }
            //std::cerr << "sumRaw=" << sum << '\n';
        }
 
        //std::cerr << "sumRaw=" << sum << '\n';
        sum = //255*pow(sum/(width*(binEnd-binStart)),1.0/power);
                255*(sum/(width*(binEnd-binStart)));
        std::cerr << "sum=" << sum << '\n';
 
 
        // TODO replace with single value!
        for (unsigned int i = 0; i < width; ++i) {
            for (unsigned int j = 0; j < height; ++j) {
                dst.at(i,j) = static_cast<T>(sum);
            }
        }
 
    };
};
 
// Computes attenuation caused by \em precipitation.
template <class T = unsigned char,class T2 = unsigned char>
class AnoRackAttenuation : public rack::AndreDetector<T,T2>
{
public:
 
    AnoRackAttenuation(const std::string & p="1.12e-4,0.62,0,0")
    : rack::AndreDetector<T,T2>("aAttenuation","Computes attenuation using cZe^p + c2Ze^p2. Alternatively, set c 'rain' or 'snow'.","c,p,c2,p2",p){
    };
 
 
    void filter(const Image<T> &src,Image<T2> &dst) const {
        const std::string & cStr = this->parameters.get("c","1.12e-4");
        double c  = this->parameters.get("c",1.12e-4);
        double p  = this->parameters.get("p",0.62);
        double c2 = this->parameters.get("c2",0);
        double p2 = this->parameters.get("p2",0);
 
        if (cStr == "rain"){
            c = 1.12E-7;
            p = 0.62;
            c2 = 0;
            p2 = 0;
        }
        else if (cStr == "snow"){
            c = 1.1E-7;
            p = 1;
            c2 = 2.1E-5;
            p2 = 0.5;
        }
        else if (cStr == "erad"){
            c = 1e-14;
            p = 0.65;
            c2 = 0.0;
            p2 = 0.0;
        }
        else if (cStr == "demo"){
            c = 1.11e-11;
            p = 0.9;
            c2 = 0.0;
            p2 = 0.0;
        }
        else if (cStr == "demo2"){
            c = 1.12e-03;
            p = 0.12;
            c2 = 0.0;
            p2 = 0.0;
        }
 
        const unsigned int srcWidth = src.getGeometry().getWidth();
        const unsigned int srcHeight = src.getGeometry().getHeight();
        const unsigned int srcChannels = src.getGeometry().getImageChannelCount();
        dst.setGeometry(srcWidth,srcHeight,srcChannels);
 
        std::cout << "Attn: << " << src.getGeometry() << "\n";
        std::cout << "Attn: << " << dst.getGeometry() << "\n";
 
        float binDepth = src.properties.get("BIN_DEPTH",500.0f);
 
        std::cerr << "Attn:" << c << ','<< p << ','<< c2 << ','<< p2 << '\n';
 
        // Main loops
 
        double dbzObs;
 
        double zAttnBin;
 
        double zAttnCumulated=0;
 
        double confidence;
        //dbzAttnCumulated;
 
        double zTrue;
 
        //double dbzTrue;
 
        for (unsigned int k = 0; k < srcChannels; ++k) {
            for (unsigned int j = 0; j < srcHeight; ++j) {
                zAttnCumulated = 0;
                for (unsigned int i = 0; i < srcWidth; ++i) {
                    dbzObs = (src.at(i,j,k)-64)*2; // FMI
 
                    zTrue = rack::dbzToZ( dbzObs ) + zAttnCumulated;
 
                    // TODO: bin width should be in the exponent?
                    zAttnBin = c*::pow(zTrue,p) * binDepth * 2.0;
                    zAttnCumulated += zAttnBin;
 
                    //confidence = (::radar::zToDbz(zAttnCumulated)+64)*4;  // TODO
                    confidence = zAttnCumulated/100.0;  // TODO
                    //confidence = 255.0 - 255.0/(1+confidence*confidence);
                    confidence = 245.0 - 245.0/(1.0+confidence);
                    //dbzTrue = ::radar::zToDbz( zTrue );
 
                    //if (j&3 == 3)                 dst.at(i,j,k) = static_cast<T2>(255.0-confidence);
                    //else
                    //dst.at(i,j,k) = static_cast<T2>(drain::radar::zToDbz(drain::radar::dbzToZ( dbzObs ) )/2 + 64);
                    //dst.at(i,j,k) = static_cast<T2>(zToDbz(zAttnCumulated));
                    dst.at(i,j,k) = static_cast<T2>(confidence);
                }
                //std::cerr << "attn = " << zAttnCumulated << "\n";
            }
        }
 
 
    }
};
 
 
template <class T = unsigned char,class T2 = unsigned char>
class AnoRack : public rack::Andre<T,T2> {
public:
 
    AnoRackShip<T,T2> ship2;
    AnoRackAttenuation<T,T2> attenuation;
    AnoRackRadomeAttenuation<T,T2> radome;
    //AnoRackRadomeAttenuation<T,T2> radome;
 
    AnoRack(){
    }
 
    virtual void addDefaultOps(){
        rack::Andre<T,T2>::addDefaultOps();
        addOperator("ship2",ship2);
        addOperator("attenuation",attenuation);
        addOperator("radome",radome);
    }
 
    virtual ~AnoRack(){};
};
 
 
}  // anorack
 
 
#endif /* ANORACK_H_ */