C++ convolve函数代码示例

void FeatureEval::blurred_intensity() {
    boost::shared_ptr<PointCloud> cloud = core_->cl_->active_;
    if(cloud == nullptr)
        return;
    int h = cloud->scan_width();
    int w = cloud->scan_height();

    boost::shared_ptr<std::vector<float>> intensity = boost::make_shared<std::vector<float>>(cloud->size());

    // Create intensity cloud
    for(uint i = 0; i < intensity->size(); i++){
        (*intensity)[i] = (*cloud)[i].intensity;
    }

    boost::shared_ptr<std::vector<float>> img = cloudToGrid(cloud->cloudToGridMap(), w*h, intensity);


    boost::shared_ptr<std::vector<float> > smooth_grad_image = convolve(img, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);
/*
    boost::shared_ptr<std::vector<float>> highfreq = img;

    for(int i = 0; i < highfreq->size(); i++){
        (*highfreq)[i] = (*highfreq)[i] - (*smooth_grad_image)[i];
    }

    drawFloats(highfreq, cloud);
*/
    drawFloats(smooth_grad_image, cloud);

}

/**
 * Hybrid window filtering, see blocks 36 and 49 of the G.728 specification.
 *
 * @param order   filter order
 * @param n       input length
 * @param non_rec number of non-recursive samples
 * @param out     filter output
 * @param hist    pointer to the input history of the filter
 * @param out     pointer to the non-recursive part of the output
 * @param out2    pointer to the recursive part of the output
 * @param window  pointer to the windowing function table
 */
static void do_hybrid_window(RA288Context *ractx,
                             int order, int n, int non_rec, float *out,
                             float *hist, float *out2, const float *window)
{
    int i;
    float buffer1[MAX_BACKWARD_FILTER_ORDER + 1];
    float buffer2[MAX_BACKWARD_FILTER_ORDER + 1];
    LOCAL_ALIGNED(32, float, work, [FFALIGN(MAX_BACKWARD_FILTER_ORDER +
                                            MAX_BACKWARD_FILTER_LEN   +
                                            MAX_BACKWARD_FILTER_NONREC, 16)]);

    av_assert2(order>=0);

    ractx->fdsp->vector_fmul(work, window, hist, FFALIGN(order + n + non_rec, 16));

    convolve(buffer1, work + order    , n      , order);
    convolve(buffer2, work + order + n, non_rec, order);

    for (i=0; i <= order; i++) {
        out2[i] = out2[i] * 0.5625 + buffer1[i];
        out [i] = out2[i]          + buffer2[i];
    }

    /* Multiply by the white noise correcting factor (WNCF). */
    *out *= 257.0 / 256.0;
}

/**
 * Hybrid window filtering. See blocks 36 and 49 of the G.728 specification.
 *
 * @param order   the order of the filter
 * @param n       the length of the input
 * @param non_rec the number of non-recursive samples
 * @param out     the filter output
 * @param in      pointer to the input of the filter
 * @param hist    pointer to the input history of the filter. It is updated by
 *                this function.
 * @param out     pointer to the non-recursive part of the output
 * @param out2    pointer to the recursive part of the output
 * @param window  pointer to the windowing function table
 */
static void do_hybrid_window(int order, int n, int non_rec, const float *in,
                             float *out, float *hist, float *out2,
                             const float *window)
{
    int i;
    float buffer1[order + 1];
    float buffer2[order + 1];
    float work[order + n + non_rec];

    /* update history */
    memmove(hist                  , hist + n, (order + non_rec)*sizeof(*hist));
    memcpy (hist + order + non_rec, in      , n                *sizeof(*hist));

    colmult(work, window, hist, order + n + non_rec);

    convolve(buffer1, work + order    , n      , order);
    convolve(buffer2, work + order + n, non_rec, order);

    for (i=0; i <= order; i++) {
        out2[i] = out2[i] * 0.5625 + buffer1[i];
        out [i] = out2[i]          + buffer2[i];
    }

    /* Multiply by the white noise correcting factor (WNCF) */
    *out *= 257./256.;
}

/**
 * Hybrid window filtering, see blocks 36 and 49 of the G.728 specification.
 *
 * @param order   filter order
 * @param n       input length
 * @param non_rec number of non-recursive samples
 * @param out     filter output
 * @param hist    pointer to the input history of the filter
 * @param out     pointer to the non-recursive part of the output
 * @param out2    pointer to the recursive part of the output
 * @param window  pointer to the windowing function table
 */
static void do_hybrid_window(int order, int n, int non_rec, float *out,
                             float *hist, float *out2, const float *window)
{
    int i;
#ifndef _MSC_VER
    float buffer1[order + 1];
    float buffer2[order + 1];
    float work[order + n + non_rec];
#else
    float *buffer1 = av_malloc_items(order + 1, float);
    float *buffer2 = av_malloc_items(order + 1, float);
    float *work = av_malloc_items(order + n + non_rec, float);
#endif

    apply_window(work, window, hist, order + n + non_rec);

    convolve(buffer1, work + order    , n      , order);
    convolve(buffer2, work + order + n, non_rec, order);

    for (i=0; i <= order; i++) {
        out2[i] = out2[i] * 0.5625 + buffer1[i];
        out [i] = out2[i]          + buffer2[i];
    }

    /* Multiply by the white noise correcting factor (WNCF). */
    *out *= 257./256.;

#ifdef _MSC_VER
	av_free(buffer1);
	av_free(buffer2);
	av_free(work);
#endif
}

/**
 * Hybrid window filtering, see blocks 36 and 49 of the G.728 specification.
 *
 * @param order   filter order
 * @param n       input length
 * @param non_rec number of non-recursive samples
 * @param out     filter output
 * @param hist    pointer to the input history of the filter
 * @param out     pointer to the non-recursive part of the output
 * @param out2    pointer to the recursive part of the output
 * @param window  pointer to the windowing function table
 */
static void do_hybrid_window(int order, int n, int non_rec, float *out,
                             float *hist, float *out2, const float *window)
{
    int i;
    #if __STDC_VERSION__ >= 199901L
    float buffer1[order + 1];
    float buffer2[order + 1];
    float work[order + n + non_rec];
    #else
    float *buffer1=(float *)alloca((order + 1)*sizeof(float));
    float *buffer2=(float *)alloca((order + 1)*sizeof(float));
    float *work=(float *)alloca((order + n + non_rec)*sizeof(float));
    #endif

    apply_window(work, window, hist, order + n + non_rec);

    convolve(buffer1, work + order    , n      , order);
    convolve(buffer2, work + order + n, non_rec, order);

    for (i=0; i <= order; i++) {
        out2[i] = out2[i] * 0.5625 + buffer1[i];
        out [i] = out2[i]          + buffer2[i];
    }

    /* Multiply by the white noise correcting factor (WNCF). */
    *out *= 257./256.;
}

  //static
  void GradientsMergeMosaic::convolveMask(weightImageType_t &rMask_p, int32 featherRadius_p)
  {
    int sideLength=1+2*featherRadius_p;
    int nElements=sideLength*sideLength;

#if 0
#if 0
    // linear filter
    LinearFilter filter(sideLength,1.0/nElements,1.0/nElements);
#else
    // box filter
    KernelFilter filter( sideLength, 1.0/nElements );
#endif
#if 1
    std::cout<<"FilterCoefficients ";
    for (const float *val=filter.Begin();val!=filter.End();++val){
      std::cout<<*val<<",";
    }
    std::cout<<std::endl;
#endif
    Convolution convolve(filter);
#else
    // separable filter
    SeparableFilter filter( sideLength,1.0/nElements);
    SeparableConvolution convolve(filter);
#endif
    convolve >>rMask_p;
  }

// TODO: THIS IS A 2D FEATURE CALCULATION! help!
// Compute the distance for each point
// Gaussian weighted average in neighbourhood
// Subtract
void FeatureEval::difference_of_gaussian_distances(){
    boost::shared_ptr<PointCloud> cloud = core_->cl_->active_;
    if(cloud == nullptr)
        return;
    int h = cloud->scan_width();
    int w = cloud->scan_height();

    // Create distance map
    boost::shared_ptr<std::vector<float>> distmap = makeDistmap(cloud);

    //distmap = interpolate(distmap, w, h, 21);

    boost::shared_ptr<std::vector<float> > smooth_grad_image = convolve(distmap, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);
    smooth_grad_image = convolve(smooth_grad_image, w, h, gaussian, 5);

    boost::shared_ptr<std::vector<float>> highfreq = distmap;

    for(uint i = 0; i < distmap->size(); i++){
        (*highfreq)[i] = (*distmap)[i] - (*smooth_grad_image)[i];
    }

    drawFloats(highfreq, cloud);
}

/*
 * Main
 */
int main(int argc, char *argv[]) {
  int i, j;
  int sobel_y[3][3] = {{-1, -2, -1},
                       {0, 0, 0},
                       {1, 2, 1}};
  int sobel_x[3][3] = {{-1, 0, 1},
                       {-2, 0, 2},
                       {-1, 0, 1}};

  Image src;
  Image G_y;
  Image G_x;
  Image out;

  if (argc < 2) {
    fprintf(stderr, "Pass the name of a PGM P5 file.\n");
    return 1;
  }

  // Read source PGM file
  if (read_PGM(argv[1], &src))
    return 1;

  // Initialize all the matrices to be of the same size as the source image.
  init_matrix(src, &G_y);
  init_matrix(src, &G_x);
  init_matrix(src, &out);

  // Convolve the Sobel filters with the source image
  convolve(sobel_y, src, &G_y);
  convolve(sobel_x, src, &G_x);

  // Calculate the gradient magnitude and put it in out.
  sobel_grad_mag(&out, G_y, G_x);

  // Normalize values to 0-255
  normalize(&G_y);
  normalize(&G_x);
  normalize(&out);

  // Write final output PGMs
  if (write_PGM("g_y.pgm", G_y))
    return 1;
  if (write_PGM("g_x.pgm", G_x))
    return 1;
  if (write_PGM("out.pgm", out))
    return 1;

  // Free the matrices
  dealloc_matrix(src.img, src.row);
  dealloc_matrix(G_y.img, 4);
  dealloc_matrix(G_x.img, 4);
  dealloc_matrix(out.img, out.row);

  return 0;
}

inline void edge_detector::edge_detect(int **a, int kernel, qreal *gradient, qreal *gangle)
{
    // Scharr kernel
    int SCx[3][3] = {{-3, 0, 3}, {-10, 0, 10}, {-3, 0, 3}};
    int SCy[3][3] = {{-3,-10,-3}, { 0, 0, 0}, { 3, 10, 3}};

    // Sobel kernel
    int Sx[3][3] = {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};
    int Sy[3][3] = {{-1,-2,-1}, { 0, 0, 0}, { 1, 2, 1}};

    // Prewitt kernel
    int Px[3][3] = {{-1, 0, 1}, {-1, 0, 1}, {-1, 0, 1}};
    int Py[3][3] = {{-1,-1,-1}, { 0, 0, 0}, { 1, 1, 1}};

    // Roberts Cross
    int Rx[3][3] = {{1, 0, 0}, {0, -1, 0}, {0, 0, 0}};
    int Ry[3][3] = {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}};

    int **Kx, **Ky;
    int k, gx, gy;

    if (kernel == 0) // Scharr kernel
    {
        k = 16;
        Kx = copy(SCx);
        Ky = copy(SCy);
    }
    else if (kernel == 1) // Sobel kernel
    {
        k = 4;
        Kx = copy(Sx);
        Ky = copy(Sy);
    }
    else if (kernel == 2) // Prewitt kernel
    {
        k = 3;
        Kx = copy(Px);
        Ky = copy(Py);
    }
    else if (kernel == 3) // Roberts Cross
    {
        k = 2;
        Kx = copy(Rx);
        Ky = copy(Ry);
    }

    gx = convolve(a, Kx, 3);
    gy = convolve(a, Ky, 3);

    deleteArray(Kx, 3);
    deleteArray(Ky, 3);

    *gradient = sqrt(gx*gx + gy*gy)/k;
    *gangle = 180/M_PI * atan2(gy, gx);
}

 static void convolve_two_segments(std::vector<point>& figure, const edge& a, const edge& b) {
   figure.clear();
   figure.push_back(point(a.first));
   figure.push_back(point(a.first));
   figure.push_back(point(a.second));
   figure.push_back(point(a.second));
   convolve(figure[0], b.second);
   convolve(figure[1], b.first);
   convolve(figure[2], b.first);
   convolve(figure[3], b.second);
 }

void image_u8_gaussian_blur(image_u8_t *im, double sigma, int ksz)
{
    assert((ksz & 1) == 1); // ksz must be odd.

    // build the kernel.
    double dk[ksz];

    // for kernel of length 5:
    // dk[0] = f(-2), dk[1] = f(-1), dk[2] = f(0), dk[3] = f(1), dk[4] = f(2)
    for (int i = 0; i < ksz; i++) {
        int x = -ksz/2 + i;
        double v = exp(-.5*sq(x / sigma));
        dk[i] = v;
    }

    // normalize
    double acc = 0;
    for (int i = 0; i < ksz; i++)
        acc += dk[i];

    for (int i = 0; i < ksz; i++)
        dk[i] /= acc;

    uint8_t k[ksz];
    for (int i = 0; i < ksz; i++)
        k[i] = dk[i]*255;

    if (0) {
        for (int i = 0; i < ksz; i++)
            printf("%d %15f %5d\n", i, dk[i], k[i]);
    }

    for (int y = 0; y < im->height; y++) {

        uint8_t x[im->stride];
        memcpy(x, &im->buf[y*im->stride], im->stride);

        convolve(x, &im->buf[y*im->stride], im->width, k, ksz);
    }

    for (int x = 0; x < im->width; x++) {
        uint8_t xb[im->height];
        uint8_t yb[im->height];

        for (int y = 0; y < im->height; y++)
            xb[y] = im->buf[y*im->stride + x];

        convolve(xb, yb, im->height, k, ksz);

        for (int y = 0; y < im->height; y++)
            im->buf[y*im->stride + x] = yb[y];
    }
}

APPROX float sobel(float w[][3]) {
    float sx;
    float sy;
    float s;

    sx = convolve(w, ky);
    sy = convolve(w, kx);

    s = sqrt(sx * sx + sy * sy);
    if (s >= (256 / sqrt(256 * 256 + 256 * 256)))
        s = 255 / sqrt(256 * 256 + 256 * 256);

    return s ;
}

/*----------------------------------------------------------------------------
 * norm_corr - Find the normalized correlation between the target vector and
 *             the filtered past excitation.
 *----------------------------------------------------------------------------
 */
static void norm_corr(
 FLOAT exc[],           /* input : excitation buffer */
 FLOAT xn[],            /* input : target vector */
 FLOAT h[],             /* input : imp response of synth and weighting flt */
 int l_subfr,           /* input : Length of frame to compute pitch */
 int t_min,             /* input : minimum value of searched range */
 int t_max,             /* input : maximum value of search range */
 FLOAT corr_norm[]      /* output: normalized correlation (correlation between
                                   target and filtered excitation divided by
                                   the square root of energy of filtered
                                    excitation) */
)
{
 int    i, j, k;
 FLOAT excf[L_SUBFR];     /* filtered past excitation */
 FLOAT  alp, s, norm;

 k = -t_min;

 /* compute the filtered excitation for the first delay t_min */

 convolve(&exc[k], h, excf, l_subfr);

 /* loop for every possible period */

 for (i = t_min; i <= t_max; i++)
 {
   /* Compute 1/sqrt(energie of excf[]) */

   alp = (F)0.01;
   for (j = 0; j < l_subfr; j++)
     alp += excf[j]*excf[j];

   norm = inv_sqrt(alp);


   /* Compute correlation between xn[] and excf[] */

   s = (F)0.0;
   for (j = 0; j < l_subfr; j++)  s += xn[j]*excf[j];


   /* Normalize correlation = correlation * (1/sqrt(energie)) */

   corr_norm[i] = s*norm;

   /* modify the filtered excitation excf[] for the next iteration */

   if (i != t_max)
   {
     k--;
     for (j = l_subfr-1; j > 0; j--)
        excf[j] = excf[j-1] + exc[k]*h[j];
     excf[0] = exc[k];
   }
 }

 return;

}

lpyramid_t *
lpyramid_create (float *image, int width, int height)
{
    lpyramid_t *pyramid;
    int i;

    pyramid = malloc (sizeof (lpyramid_t));
    if (pyramid == NULL) {
	fprintf (stderr, "Out of memory.\n");
	exit (1);
    }
    pyramid->width = width;
    pyramid->height = height;

    /* Make the Laplacian pyramid by successively
     * copying the earlier levels and blurring them */
    for (i=0; i<MAX_PYR_LEVELS; i++) {
	pyramid->levels[i] = malloc (width * height * sizeof (float));
	if (pyramid->levels[i] == NULL) {
	    fprintf (stderr, "Out of memory.\n");
	    exit (1);
	}
	if (i == 0) {
	    memcpy (pyramid->levels[i], image, width * height * sizeof (float));
	} else {
	    convolve(pyramid, pyramid->levels[i], pyramid->levels[i - 1]);
	}
    }

    return pyramid;
}

void CONVOLVE1::process(const float *buf, const int len)
{
	DPRINT1("CONVOLVE1::process (len=%d)\n", len);

	// NOTE: <len> will always be equal to _impframes
	if (_winosc)
		for (int i = 0; i < len; i++)
			_fftbuf[i] = buf[i] * _winosc->next(i);
	else
		for (int i = 0; i < len; i++)
			_fftbuf[i] = buf[i];
	for (int i = len; i < _fftlen; i++)
		_fftbuf[i] = 0.0f;

	convolve();

	for (int i = 0, j = len; i < len; i++, j++) {
		_ovadd[i] += _fftbuf[i];
		_wet[_outWriteIndex] = _ovadd[i];
		_dry[_outWriteIndex] = buf[i];
		incrementOutWriteIndex();
		_ovadd[i] = _fftbuf[j];
	}

//	for (int i = 0; i < _fftlen; i++) printf("[%d] = %f\n", i, _fftbuf[i]);
}