C++ Mat_::mul方法代码示例

//===========================================================================
void Multi_SVR_patch_expert::Response(const Mat_<float> &area_of_interest, Mat_<double> &response)
{
	
	int response_height = area_of_interest.rows - height + 1;
	int response_width = area_of_interest.cols - width + 1;

	if(response.rows != response_height || response.cols != response_width)
	{
		response.create(response_height, response_width);
	}

	// For the purposes of the experiment only use the response of normal intensity, for fair comparison

	if(svr_patch_experts.size() == 1)
	{
		svr_patch_experts[0].Response(area_of_interest, response);		
	}
	else
	{
		// responses from multiple patch experts these can be gradients, LBPs etc.
		response.setTo(1.0);
		
		Mat_<double> modality_resp(response_height, response_width);

		for(size_t i = 0; i < svr_patch_experts.size(); i++)
		{			
			svr_patch_experts[i].Response(area_of_interest, modality_resp);			
			response = response.mul(modality_resp);	
		}	
		
	}

}

inline float calculate_var(const Mat_<float>& v1){
    if (v1.rows==0)
    {
		return 0;
    }
	float mean_1 = mean(v1)[0];
    float mean_2 = mean(v1.mul(v1))[0];
    return mean_2 - mean_1*mean_1;
    
}

void iterative_computation(Mat_<float>& u, Mat_<float>& v, const Mat_<float>& Ix, 
			   const Mat_<float>& Iy, const Mat_<float>& It) {
	if (METHOD == 0) { // Jacobi Methd
	    Mat_<float> f =  (Mat_<float>(3,3) << 1.0/12.0, 1.0/6.0, 1.0/12.0, 1.0/6.0, 0.0, 1.0/6.0, 
					  1.0/12.0, 1.0/6.0, 1.0/12.0);
	    Mat_<float> avg_u, avg_v;
	    filter2D(u, avg_u, -1 , f, Point(-1, -1), 0, BORDER_DEFAULT );
	    filter2D(v, avg_v, -1 , f, Point(-1, -1), 0, BORDER_DEFAULT );
	    Mat_<float> d1 = Ix.mul(avg_u) + Iy.mul(avg_v) + It;
	    Mat_<float> d2 = Mat::ones(u.size(), CV_32F) * ALPHA * ALPHA + Ix.mul(Ix) + Iy.mul(Iy);
	    Mat_<float> r = d1.mul(1 / d2);
	    u = avg_u - Ix.mul(r);
	    v = avg_v - Iy.mul(r);
	}
	else if (METHOD == 1) { // Gauss-Seidel method 
	    for (int i = 1; i < u.rows - 1; i ++) {
		for (int j = 1; j < u.cols - 1; j ++) {
			float avg_u = 1.f / 6.f * (u.at<float>(i - 1, j) + u.at<float>(i, j + 1) + 
					       u.at<float>(i + 1, j) + u.at<float>(i, j - 1)) +
				      1.f / 12.f * (u.at<float>(i - 1, j -1) + u.at<float>(i - 1, j + 1) + 
						u.at<float>(i + 1, j - 1) + u.at<float>(i + 1, j + 1));
			float avg_v = 1.f / 6.f * (v.at<float>(i - 1, j) + v.at<float>(i, j + 1) + 
					       v.at<float>(i + 1, j) + v.at<float>(i, j - 1)) +
				      1.f / 12.f * (v.at<float>(i - 1, j -1) + v.at<float>(i - 1, j + 1) + 
						v.at<float>(i + 1, j - 1) + v.at<float>(i + 1, j + 1));
			float ix = Ix.at<float>(i, j);
			float iy = Iy.at<float>(i, j);
			float it = It.at<float>(i, j);
			float r = (ix * avg_u + iy * avg_v + it) / 
				  (ALPHA * ALPHA + ix * ix + iy * iy);
			
			u.at<float>(i, j) = avg_u - ix * r;
			v.at<float>(i, j) = avg_v - iy * r;
		}
	    }
	} 

	else if (METHOD == 2) { // Successive overrelaxation method 
	    for (int i = 1; i < u.rows - 1; i ++) {
		for (int j = 1; j < u.cols - 1; j ++) {
			float avg_u = 1.f / 6.f * (u.at<float>(i - 1, j) + u.at<float>(i, j + 1) + 
					       u.at<float>(i + 1, j) + u.at<float>(i, j - 1)) +
				      1.f / 12.f * (u.at<float>(i - 1, j -1) + u.at<float>(i - 1, j + 1) + 
						u.at<float>(i + 1, j - 1) + u.at<float>(i + 1, j + 1));
			float avg_v = 1.f / 6.f * (v.at<float>(i - 1, j) + v.at<float>(i, j + 1) + 
					       v.at<float>(i + 1, j) + v.at<float>(i, j - 1)) +
				      1.f / 12.f * (v.at<float>(i - 1, j -1) + v.at<float>(i - 1, j + 1) + 
						v.at<float>(i + 1, j - 1) + v.at<float>(i + 1, j + 1));
			float ix = Ix.at<float>(i, j);
			float iy = Iy.at<float>(i, j);
			float it = It.at<float>(i, j);
			float r = (ix * avg_u + iy * avg_v + it) / 
				  (ALPHA * ALPHA + ix * ix + iy * iy);
			
			u.at<float>(i, j) = (1 - WEIGHT) * u.at<float>(i, j) + WEIGHT * (avg_u - ix * r);
			v.at<float>(i, j) = (1 - WEIGHT) * v.at<float>(i, j) + WEIGHT * (avg_v - iy * r);
		}
	    }
	}			   
}  

/**
 * @brief       Window Function is just an element-wise product, a type of weighted staff. 
 *              Absolutely not a convolution !!!
 * @param       iImg    Input   -- the original signal
 * @param       oImg    Output  -- the transformed signal
 */
void VO_STFT::VO_ForwardTransform(  const Mat_<float>& iImg,
                                    Mat_<float>& oImg) const
{
    if( (iImg.cols != this->m_VOWindowFunc->m_MatWindowedKernel.cols) ||
        (iImg.rows != this->m_VOWindowFunc->m_MatWindowedKernel.rows) )
    {
        cerr << "STFT: signal should have the same size as Window kernel" << endl;
        exit(1);
    }
    // Explained by JIA Pei. for STFT, it's element-wise weighting, rather than convolution!!
    oImg = iImg.mul(this->m_VOWindowFunc->m_MatWindowedKernel);
    cv::dft(oImg, oImg);
}

inline double calculate_var(const Mat_<double>& v1){
	double mean_1 = mean(v1)[0];
	double mean_2 = mean(v1.mul(v1))[0];
	return mean_2 - mean_1*mean_1;
	
}