2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 свертка 3D-массивов на С++
Сообщение02.09.2015, 22:02 
Экс-модератор
Аватара пользователя


23/12/05
11256
Не уверен, что этот раздел будет самым подходящим, но попробуем тут.
Стоит задача перепилить кое-какой матлабовский код на С++. Среди прочего в коде столкнулся с функцией convn(X, Y, shape);
В моем случае X,Y - трехмерные не очень большие массивы (могут быть разные, но то, что смотрел для одного из случаев было что-то типа 3х3х27 и 1х1х3, но могут быть и другие), shape= 'same'.

Буду признателен, если кто-то поможет разобраться с алгоритмом - что конкретно делает эта функция и с реализацией на плюсах.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение02.09.2015, 22:49 
Аватара пользователя


31/10/08
1242
Вот тут наглядно показана свёртка для одномерного случая и для двухмерного, для 3-х мерного по аналогии. Для N-мерного по индукции.
http://www.songho.ca/dsp/convolution/convolution.html

Ядро перемещается в пространстве исходных данных. Исходные данные в окрестности ядра перемножаются и суммируются. Результат записывается в результирующий массив.

shape - задаёт режим обработки границ.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 11:03 
Экс-модератор
Аватара пользователя


23/12/05
11256
Pavia в сообщении #1050029 писал(а):
Вот тут наглядно показана свёртка для одномерного случая и для двухмерного, для 3-х мерного по аналогии. Для N-мерного по индукции. http://www.songho.ca/dsp/convolution/convolution.html
У меня сайт не открывается ((

Pavia в сообщении #1050029 писал(а):
Ядро перемещается в пространстве исходных данных. Исходные данные в окрестности ядра перемножаются и суммируются.

В целом-то, я примерно представляю, но, видимо, пространственного воображения не хватает, чтобы правильно расписать индексы массивов, что с чем перемножается и суммируется...

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 12:32 
Заслуженный участник
Аватара пользователя


30/01/06
09/12/19
71257
http://www.mathworks.com/help/matlab/ref/conv2.html
хотя бы открывается?

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 14:11 
Экс-модератор
Аватара пользователя


23/12/05
11256
Munin в сообщении #1050118 писал(а):
http://www.mathworks.com/help/matlab/ref/conv2.html
хотя бы открывается?

Да, это открывается )

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 16:16 
Заслуженный участник
Аватара пользователя


30/01/06
09/12/19
71257
Тут, я так понимаю, нам предлагается самостоятельно обобщить формулы
$$c(i)=\sum\limits_j a(j)b(i-j+1)\qquad\text{и}\qquad c(n_1,n_2)=\sum\limits_{k_1=-\infty}^\infty\sum\limits_{k_2=-\infty}^\infty a(k_1,k_2)b(n_1-k_1,n_2-k_2)$$ на случай $N$ измерений. Сделать это нетрудно даже без пространственного воображения, надо только аккуратно писать формулы. Единственную трудность я вижу в точном понимании того, что именно подразумевается под "центральной частью" результата, в случае, когда соответствующие размерности массивов чётные или нечётные, и поэтому не делятся точно на 2.

Если это лучшая доступная документация, то дальше два пути: либо добывать исходники, либо выяснять поведение эталонной системы экспериментами. Нормально, для портирования часто так бывает :-)

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 17:54 
Экс-модератор
Аватара пользователя


23/12/05
11256
Если я правильно представил и записал, то для массивов $X$ размера $a_1\times a_2\times a_3$ и $Y$ размера $b_1\times b_2\times b_3$ свертка
$Z$ размера $c_1\times c_2\times c_3$ (где $c_i=a_i+b_i-1$) будет такой:

$Z(k_1,k_2,k_3)
=\sum\limits_{j_1=\max(1,k_1+1-b_1)}^{\min(k_1,a_1)}\sum\limits_{j_2=\max(1,k_2+1-b_2)}^{\min(k_2,a_2)}\sum\limits_{j_3=\max(1,k_3+1-b_3)}^{\min(k_3,a_3)}X(j_1,j_2,j_3)Y(k_1+1-j_1,k_2+1-j_2,k_3+1-j_3)$


А потом выбрать центральную часть размера $a_1\times a_2\times a_3$
В случае, если четность $a_i$, $c_i$ не совпадет, просто выбрать какое-то предпочтительное направление от центра и округлить либо в меньшую, либо в большую сторону положение центра.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 18:35 
Заслуженный участник
Аватара пользователя


30/01/06
09/12/19
71257
Чё-то трудно мне с налёту и без комментариев в этих минимумах и максимумах разобраться.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 18:46 
Экс-модератор
Аватара пользователя


23/12/05
11256
Писал по аналогии с описанием для одномерного случая отсюда: http://www.exponenta.ru/soft/matlab/potemkin/book2/chapter8/conv.asp

Сейчас заимплементить уже не успею - надо бежать, а завтра сделаю, сравню результаты с матлабовскими... Попутно, можно и 1D, 2D сделать - их проще реализовать и проверить.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 18:55 
Заслуженный участник
Аватара пользователя


30/01/06
09/12/19
71257
photon в сообщении #1050214 писал(а):
Писал по аналогии с описанием для одномерного случая отсюда: http://www.exponenta.ru/soft/matlab/potemkin/book2/chapter8/conv.asp

Ну, я бы так сразу не списывал из чёрт-знает-какого источника.

photon в сообщении #1050214 писал(а):
Сейчас заимплементить уже не успею - надо бежать, а завтра сделаю, сравню результаты с матлабовскими...

Тестирование продумайте.

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 19:00 
Заслуженный участник
Аватара пользователя


20/01/06
1030
Я бы основывался на этой штуке http://rosettacode.org/wiki/Image_convolution,
а вот исходник из матлаба (октава) http://code.google.com/p/addi/source/browse/Addi/assets/toolbox/polynomial/convn.m?r=87

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 19:20 
Заслуженный участник
Аватара пользователя


04/09/14
3859
ФТИ им. Иоффе СПб
Я, конечно, не в своем огороде копаюсь, но... Команда convn (исходник лежит в \MATLAB701\toolbox\matlab\datafun\convn.m) использует convnc (лежит в \MATLAB701\toolbox\matlab\datafun\convnc.cpp), которая написана на С, поэтому, IMHO, исходная проблема переписывания отсутствует. (Если что не так - заранее извиняюсь).

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 19:32 
Заслуженный участник
Аватара пользователя


30/01/06
09/12/19
71257
Freude в сообщении #1050224 писал(а):
а вот исходник из матлаба (октава) http://code.google.com/p/addi/source/browse/Addi/assets/toolbox/polynomial/convn.m?r=87

Вот что там лежит:
:-)

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 20:01 
Экс-модератор
Аватара пользователя


23/12/05
11256
Freude в сообщении #1050224 писал(а):
а вот исходник из матлаба (октава) http://code.google.com/p/addi/source/br ... nvn.m?r=87


Это я видел.
amon в сообщении #1050232 писал(а):
(лежит в \MATLAB701\toolbox\matlab\datafun\convnc.cpp), которая написана на С

А это, спасибо, проверю, когда теперь доберусь до матлаба. Если действительно там есть исходники, то, конечно, это будет идеальный вариант... А если нет, то попрошу у вас ))

 Профиль  
                  
 
 Re: свертка 3D-массивов на С++
Сообщение03.09.2015, 20:26 
Заслуженный участник
Аватара пользователя


04/09/14
3859
ФТИ им. Иоффе СПб
photon в сообщении #1050250 писал(а):
А если нет, то попрошу у вас

Не откладывая и на всякий случай (меня может не быть завтра).
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
/*

   CONVNC.CPP   .MEX file
   Implements full N-D convolution.
   Inputs must be real, numeric, and float.

   Only one syntax is supported:
     C = CONVNC(A,B)

   Copyright 1984-2004 The MathWorks, Inc.
*/


#include "mex.h"

static char rcsid[] = "$Revision $";

/* Input and output arguments */
#define A (prhs[0])
#define B (prhs[1])
#define C (plhs[0])

/* Increment subscripts */
/* Increment subscript vector by one element */
/* in the direction of linear indexing. */
/* The subscript vector is assumed to be zero-based. */

inline void INCREMENT_SUBSCRIPTS(int *subs, const int*size, int ndims) {
  subs[0] += 1;
  for (int p = 0; p < (ndims-1); p++) {
    if (subs[p] > (size[p]-1)) {
      subs[p] = 0;
      subs[p+1] += 1;
    }
  }
}

/* Convert a zero-based subscript vector to a linear index. */

inline void SUBSCRIPTS_TO_LINEAR(const int *subs, int ndims, const int *pdims,
                                 int *index)
{
  int i = ndims;
  int factor = 1;
                 
  *index = 0;
             
  while (i--) {
    *index += *subs++ * factor;
    factor *= *pdims++;
  }
}

template <class Tout, class Tin1, class Tin2>
void convolve(Tout *c, const Tin1 *a, const Tin2 *b,
                     const int *sizeA, const int *sizeB, const int *sizeC,
                     int ndimsA, int ndimsB, int ndimsC)
{
  /* Input sanity checking */
  /* Any of the following error messages indicates that */
  /* something went seriously wrong in mexFunction().   */
  /* I would use utAssert() if the API had it.  -sle    */
  if (a==NULL || b==NULL || c==NULL) {
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError1",
    "Internal consistency error in convolve().");
  }
  if (sizeA==NULL || sizeB==NULL || sizeC==NULL) {
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError1",
    "Internal consistency error in convolve().");
  }
  if ((ndimsA != ndimsB) || (ndimsA != ndimsC)) {
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError1",
    "Internal consistency error in convolve().");
  }

  /* Compute the number of elements in inputs a and b */
  int lengthA = 1;     /* length of input A */
  for (int p = 0; p < ndimsA; p++) {
    lengthA *= sizeA[p];
  }
  int lengthB = 1;     /* length of input B */
  for (int p = 0; p < ndimsB; p++) {
    lengthB *= sizeB[p];
  }

  /* Initialize subscript vectors */
  int *subsA = (int *) mxCalloc(ndimsA, sizeof(int)); /* subscript vector for A */
  int *subsB = (int *) mxCalloc(ndimsB, sizeof(int)); /* subscript vector for B */
  int *subsC = (int *) mxCalloc(ndimsC, sizeof(int)); /* subscript vector for C */
  if (subsA==NULL || subsB==NULL || subsC==NULL) {
    if (subsA != NULL) {
      mxFree(subsA);
    }
    if (subsB != NULL) {
      mxFree(subsB);
    }
    if (subsC != NULL) {
      mxFree(subsC);
    }
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
    "Memory allocation failure.");
  }

  /* Initialize subscript array for a to be a(-1,0,0,...,0) */
  /* This is ok since subscripts will be incremented once before use. */
  subsA[0] = -1;
  for (int p = 1; p < ndimsA; p++) {
    subsA[p] = 0;
  }

  /* Core computation loops */
  int linearIndexC=0;  /* linear index into output array */
  for (int p = 0; p < lengthA; p++) {
    /* Increment subscript vector for A */
    INCREMENT_SUBSCRIPTS(subsA, sizeA, ndimsA);

    /* Initialize subscript array for a to be b(-1,0,0,...,0) */
    /* This is ok since subscripts will be incremented once before use. */
    subsB[0] = -1;
    for (int r = 1; r < ndimsA; r++) {
      subsB[r] = 0;
    }
    for (int q = 0; q < lengthB; q++) {
      /* Increment subscript vector for B */
      INCREMENT_SUBSCRIPTS(subsB, sizeB, ndimsB);
     
      /* Where should the next partial product go in the output array? */
      /* Answer: subsC = subsA + subsB */
      for (int r = 0; r < ndimsA; r++) {
        subsC[r] = subsA[r] + subsB[r];
      }

      /* But we need the answer as a linear index rather than a */
      /* subscript array */
      SUBSCRIPTS_TO_LINEAR(subsC, ndimsC, sizeC, &linearIndexC);
 
      /* Accumulate partial product */
      c[linearIndexC] += a[p] * b[q];
    }
  }

  /* Clean up and go home */
  mxFree(subsA);
  mxFree(subsB);
  mxFree(subsC);
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  const int NumInputs = 2;

  /* Sanity check input */
  if (nrhs != NumInputs) mexErrMsgIdAndTxt("MATLAB:convnc:WrongNumInputs",
    "Two input arguments required.");

  bool singleA = mxIsSingle(A);
  bool singleB = mxIsSingle(B);
 
  if ((!mxIsDouble(A) && !singleA) || mxIsComplex(A) ||
      (!mxIsDouble(B) && !singleB) || mxIsComplex(B)) {
     mexErrMsgIdAndTxt("MATLAB:convnc:UnsupportedDataType",
    "Inputs must be real and float (double or single).");
  }
 
  mxClassID outputClass = (singleA || singleB ? mxSINGLE_CLASS : mxDOUBLE_CLASS);

  /* If either A or B is empty, return an empty matrix and go home. */
  if (mxIsEmpty(A) || mxIsEmpty(B)) {
    C = mxCreateNumericMatrix(0,0, outputClass, mxREAL);
    return;
  }

  int ndimsA = mxGetNumberOfDimensions(A); /* Dimensionality of first input */
  int ndimsB = mxGetNumberOfDimensions(B); /* Dimensionality of second input */
  const int *sizeA = mxGetDimensions(A);   /* Size of first input */
  const int *sizeB = mxGetDimensions(B);   /* Size of second input */

  /* Make dimensionality of A and B conform */
  int *adjustedSizeA=NULL;  /* Adjusted size of first input */
  int *adjustedSizeB=NULL;  /* Adjusted size of second input */
  int adjustedNDimsA=0;     /* Adjusted dimensionality of first input */
  int adjustedNDimsB=0;     /* Adjusted dimensionality of second input */
  if (ndimsA != ndimsB) {
    if (ndimsA > ndimsB) {
      adjustedNDimsA = ndimsA;
      adjustedNDimsB = ndimsA;
      adjustedSizeB = (int *) mxCalloc(adjustedNDimsB, sizeof(int));
      adjustedSizeA = (int *) mxCalloc(adjustedNDimsA, sizeof(int));
      if (adjustedSizeB == NULL || adjustedSizeA == NULL) {
        if (adjustedSizeB != NULL) mxFree((void *) adjustedSizeB);
        if (adjustedSizeA != NULL) mxFree((void *) adjustedSizeA);
        mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
                          "Memory allocation failure.");
      }
      for (int p = 0; p < ndimsB; p++) {
        adjustedSizeB[p] = sizeB[p];
        adjustedSizeA[p] = sizeA[p];
      }
      for (int p = ndimsB; p < adjustedNDimsB; p++) {
        adjustedSizeB[p] = 1;
        adjustedSizeA[p] = sizeA[p];
      }
    } else {
      adjustedNDimsA = ndimsB;
      adjustedNDimsB = ndimsB;
      adjustedSizeA = (int *) mxCalloc(adjustedNDimsA, sizeof(int));
      adjustedSizeB = (int *) mxCalloc(adjustedNDimsB, sizeof(int));
      if (adjustedSizeA == NULL || adjustedSizeB == NULL) {
        if (adjustedSizeA != NULL) mxFree((void *) adjustedSizeA);
        if (adjustedSizeB != NULL) mxFree((void *) adjustedSizeB);
        mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
                          "Memory allocation failure.");
      }
      for (int p = 0; p < ndimsA; p++) {
        adjustedSizeA[p] = sizeA[p];
        adjustedSizeB[p] = sizeB[p];
      }
      for (int p = ndimsA; p < adjustedNDimsA; p++) {
        adjustedSizeA[p] = 1;
        adjustedSizeB[p] = sizeB[p];
      }
    }
  } else {
    adjustedNDimsA = ndimsA;
    adjustedNDimsB = ndimsB;
    adjustedSizeA = (int *) mxCalloc(adjustedNDimsA, sizeof(int));
    adjustedSizeB = (int *) mxCalloc(adjustedNDimsB, sizeof(int));
    if (adjustedSizeA == NULL || adjustedSizeB == NULL) {
      if (adjustedSizeA != NULL) mxFree((void *) adjustedSizeA);
      if (adjustedSizeB != NULL) mxFree((void *) adjustedSizeB);
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
                      "Memory allocation failure.");
    }
    for (int p = 0; p < adjustedNDimsA; p++) {
      adjustedSizeA[p] = sizeA[p];
      adjustedSizeB[p] = sizeB[p];
    }
  }

  /* Initialize output */
  int *sizeC = (int *) mxCalloc(adjustedNDimsA, sizeof(int));
  if (sizeC == NULL) {
    mxFree((void *) adjustedSizeA);
    mxFree((void *) adjustedSizeB);
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
                      "Memory allocation failure.");
  }
  for (int p = 0; p < adjustedNDimsA; p++) {
    sizeC[p] = adjustedSizeA[p] + adjustedSizeB[p] - 1;
  }
  C = mxCreateNumericArray(adjustedNDimsA, sizeC, outputClass, mxREAL);
  if (C == NULL) {
    mxFree((void *) adjustedSizeA);
    mxFree((void *) adjustedSizeB);
    mxFree((void *) sizeC);
    mexErrMsgIdAndTxt("MATLAB:convnc:InternalError2",
                      "Memory allocation failure.");
  }

  if (singleA && singleB) {  

    real32_T *sA = (real32_T *) mxGetData(A);
    real32_T *sB = (real32_T *) mxGetData(B);
    real32_T *sC = (real32_T *) mxGetData(C);
    convolve(sC, sA, sB,
           adjustedSizeA, adjustedSizeB, sizeC,
           adjustedNDimsA, adjustedNDimsB, adjustedNDimsB);

  } else  if (singleA && !singleB) {

    real32_T *sA = (real32_T *) mxGetData(A);
    double   *dB = (double   *) mxGetData(B);
    real32_T *sC = (real32_T *) mxGetData(C);
    convolve(sC, sA, dB,
             adjustedSizeA, adjustedSizeB, sizeC,
             adjustedNDimsA, adjustedNDimsB, adjustedNDimsB);

  } else if (!singleA && singleB) {

    double   *dA = (double   *) mxGetData(A);
    real32_T *sB = (real32_T *) mxGetData(B);
    real32_T *sC = (real32_T *) mxGetData(C);
    convolve(sC, dA, sB,
           adjustedSizeA, adjustedSizeB, sizeC,
           adjustedNDimsA, adjustedNDimsB, adjustedNDimsB);

  } else {    /* both inputs are double */

    double *dA = (double *) mxGetData(A);
    double *dB = (double *) mxGetData(B);  
    double *dC = (double *) mxGetData(C);
    convolve(dC, dA, dB,
           adjustedSizeA, adjustedSizeB, sizeC,
           adjustedNDimsA, adjustedNDimsB, adjustedNDimsB);

  }
}
 

А это - код conv
код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
function C = convn(A,B,shape)
%CONVN  N-dimensional convolution.
%   C = CONVN(A, B) performs the N-dimensional convolution of
%   matrices A and B.
%   C = CONVN(A, B, 'shape') controls the size of the answer C:
%     'full'   - (default) returns the full N-D convolution
%     'same'   - returns the central part of the convolution that
%                is the same size as A.
%     'valid'  - returns only the part of the result that can be
%                computed without assuming zero-padded arrays.  The
%                size of the result is max(size(A)-size(B)+1,0).
%
%   Class support for inputs A,B:
%      float: double, single
%
%   See also CONV, CONV2.

%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 1.10.4.3 $  $Date: 2004/03/09 16:16:09 $

error(nargchk(2,3,nargin));

if (nargin < 3)
  shape = 'full';
end

if (issparse(A) || issparse(B))
  error('MATLAB:convn:SparseInput', 'Input arguments must be full.')
end

if ~isfloat(A)
    A = double(A);
end
if ~isfloat(B)
    B = double(B);
end

Aisreal = isreal(A);
Bisreal = isreal(B);

if (Aisreal && Bisreal)
  C = convnc(A,B);
 
elseif (Aisreal && ~Bisreal)
  C = convnc(A,real(B)) + j*convnc(A,imag(B));
 
elseif (~Aisreal && Bisreal)
  C = convnc(real(A),B) + j*convnc(imag(A),B);
 
else
  Ar = real(A);
  Ai = imag(A);
  Br = real(B);
  Bi = imag(B);
  C = convnc(Ar,Br) - convnc(Ai,Bi) + j*(convnc(Ai,Br) + convnc(Ar,Bi));
 
end

if (strcmp(shape,'same'))
  sizeA = [size(A) ones(1,ndims(C)-ndims(A))];
  sizeB = [size(B) ones(1,ndims(C)-ndims(B))];
  flippedKernelCenter = ceil((1 + sizeB)/2);
  subs = cell(1,ndims(C));
  for p = 1:length(subs)
    subs{p} = (1:sizeA(p)) + flippedKernelCenter(p) - 1;
  end
  C = C(subs{:});
 
elseif (strcmp(shape,'valid'))
  sizeB = [size(B) ones(1,ndims(C)-ndims(B))];
  outSize = max([size(A) ones(1,ndims(C)-ndims(A))] - sizeB + 1, 0);
  subs = cell(1,ndims(C));
  for p = 1:length(subs)
    subs{p} = (1:outSize(p)) + sizeB(p) - 1;
  end
  C = C(subs{:});
 
end
 

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 22 ]  На страницу 1, 2  След.

Модераторы: Toucan, maxal, PAV, Karan, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group