cv::Mat MalikPerona_quadr(cv::Mat _src, float k)
        {
            cv::Mat dxM ;
            cv::Mat dyM;
            Sobel(_src, dxM, CV_32F, 1, 0, 3);
            Sobel(_src, dyM, CV_32F, 0, 1, 3);

            cv::Mat cM = cv::Mat::zeros(cv::Size(_src.cols, _src.rows), CV_32FC1);
            for(int i = 0; i < _src.cols; i++)
                for(int j = 0; j < _src.rows; j++)
                    cM.at<float>(j,i) = 1.0f / (1.0f + (std::pow(dxM.at<float>(j,i), 2) + std::pow(dyM.at<float>(j,i), 2))/k*k );
            return cM;
        }

        cv::Mat MalikPerona_exp(cv::Mat _src, float k)
        {
            cv::Mat dxM ;
            cv::Mat dyM;
            Sobel(_src, dxM, CV_32F, 1, 0, 3);
            Sobel(_src, dyM, CV_32F, 0, 1, 3);

            cv::Mat cM = cv::Mat::zeros(cv::Size(_src.cols, _src.rows), CV_32FC1);
            for(int i = 0; i < _src.cols; i++)
                for(int j = 0; j < _src.rows; j++){
                    cM.at<float>(j,i) = std::exp(-((std::pow(dxM.at<float>(j,i), 2) + std::pow(dyM.at<float>(j,i), 2))/0.05f*0.05f));
                }
            return cM;

        }

        cv::Mat AnisotropicDiffusion(cv::Mat _src, int type, float k , float dt, float lambda )
        {
            _src.convertTo(_src, CV_32FC1, 1.0/255.0);
            cv::Mat cM;
            if(type == EXP)
                cM = MalikPerona_exp(_src, k);
            else if(type == QUADR)
                cM = MalikPerona_quadr(_src, k);
            else
                cM = prop(_src, k);

            cv::Mat dst = cv::Mat::zeros(cv::Size(_src.cols, _src.rows), CV_32FC1);

            for(int i = 1; i < _src.cols - 1; i++)
                for(int j = 1; j < _src.rows - 1; j++){

                    float v = ((cM.at<float>(j,i-1) + cM.at<float>(j,i))*(_src.at<float>(j,i-1) - _src.at<float>(j,i)) +
                                (cM.at<float>(j,i+1) + cM.at<float>(j,i))*(_src.at<float>(j,i+1) - _src.at<float>(j,i)) +
                                (cM.at<float>(j-1,i) + cM.at<float>(j,i))*(_src.at<float>(j-1,i) - _src.at<float>(j,i)) +
                                (cM.at<float>(j+1,i) + cM.at<float>(j,i))*(_src.at<float>(j+1,i) - _src.at<float>(j,i)));
                    dst.at<float>(j,i) = _src.at<float>(j,i) + dt*lambda*v;
                }
            _src.convertTo(_src, CV_8UC1, 255.0/1.0);
            dst.convertTo(dst, CV_8UC1, 255.0/1.0);

            return dst;
        }