Проблемы производительности алгоритма взаимной корреляции C # Bitmap - PullRequest
2 голосов
/ 23 марта 2012

Это мой первый вопрос о переполнении стека, так что я нашел все ответы на свои проблемы в кратчайшие сроки. Большое спасибо! Обычно я в основном программирую на ПЛК, мои знания о мире ПК весьма ограничены, и я впервые использую C #.

Итак, мне довелось попытаться взаимно коррелировать две области пикселей в двух растровых изображениях, согласно статье здесь: http://users.ox.ac.uk/~atdgroup/publications/Rankov,%20V.,%20Proceedings%20of%20Spie,%20Vol.%205701,2005.pdf

[EDIT] Цель состоит в том, чтобы найти точное местоположение совпадения, чтобы выполнить сшивание двух изображений. Я также вынул часть прокомментированного кода, чтобы сделать обзор лучше (я открою еще один вопрос для части скользящего среднего). [/ EDIT]

Мои проблемы - правильная реализация скользящей средней и общие настройки производительности, где, я надеюсь, вы, ребята, сможете мне помочь.

Растровые изображения имеют фиксированное перекрытие во всех известных мне направлениях (10%), поэтому я могу сохранять области поиска (называемые составной областью в исходном коде ниже) довольно небольшими, но не настолько маленькими, как кажется. Я также предполагаю, что они имеют одинаковый размер и формат пикселей. Однако производительность моего алгоритма меня не удовлетворяет. У меня есть чувство (в основном из-за того, что мне не хватает «глубоких» знаний и опыта), что есть много возможностей для совершенствования.

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

  • Расчет значений пикселей в отдельном методе (в основном для удобства чтения, быстро отбрасывается)
  • четыре вложенных в петли

Вот некоторые временные измерения выпуска (Core Duo 2,4 ГГц, 4 ГБ) для двух 950px * 950px, 24RGB растровых изображений. Область поиска (область составного изображения) составляла 70px * 800px, область выборки 8px * 400px.

  • отдельная средняя функция: 5519мс
  • средняя встроенная функция: 5350 мс (только?)
  • [EDIT] изменения, предложенные Yaur: 700 мс! [/ EDIT]

Как правило, использование небольших выборок и областей поиска (4x40 и 30x100) дает довольно быстрое время, варьирующееся от нескольких мс до нескольких мс. К сожалению, чтобы быть в безопасности, чтобы найти совпадения, я должен использовать большие площади. Прежде чем перейти к субсэмплингу и т. Д., Я бы хотел быть уверенным, что мой нынешний алгоритм не полностью вне мира.

Есть какие-нибудь хитрости / уловки или общие улучшения, о которых вы можете подумать? Любая подсказка будет с благодарностью оценена.

[EDIT] Метод корреляции (значительно улучшен):

private unsafe void CrossCorrelate(ref float CCCoefficient, ref Point SampleMatchLocation)
{
    float res = 0;
    float tmpRes = 0;

    // get bit data of sample area
    BitmapData bmdSample = m_bmpSampleRaw.LockBits(m_rectSampleArea, ImageLockMode.ReadOnly, m_bmpSampleRaw.PixelFormat);
    byte* pSample = (byte*)(void*)bmdSample.Scan0;

    // calculate sample average and coefficient 1 (stays same for all iterations)
    int SampleAvg = GetAverage(bmdSample, 0, bmdSample.Width, 0, bmdSample.Height);
    float CN1     = GetCN1(bmdSample, SampleAvg);

    int CompAvg         = 0;
    BitmapData bmdComp  = null;
    Rectangle compRect;

    int SearchHeightLimit   = m_rectSearchArea.Height - m_rectSampleArea.Height;
    int SearchWidthLimit    = m_rectSearchArea.Width - m_rectSampleArea.Width;
    int SearchLocX          = m_rectSearchArea.X;
    int SearchLocY          = m_rectSearchArea.Y;
    int SampleHeight        = m_rectSampleArea.Height;
    int SampleWidth         = m_rectSampleArea.Width;

    int a = 0; // used to calculate power of 2 without using Math.Pow

    // iterate through search area, 
    // in case of equal sizes make sure it iterates at least once
    if (SearchHeightLimit == 0) SearchHeightLimit++;
    if (SearchWidthLimit == 0) SearchWidthLimit++;

    for (int i = 0; i < SearchHeightLimit; i++)
    {
        for (int j = 0; j < SearchWidthLimit; j++)
        {
            int CN0Sum = 0;
            int CN2Sum = 0;

            // create composite pixel data at current search location
            compRect    = new Rectangle(SearchLocX + j, SearchLocY + i, SampleWidth, SampleHeight);
            bmdComp     = m_bmpCompositeRaw.LockBits(compRect, ImageLockMode.ReadOnly, m_bmpCompositeRaw.PixelFormat);
            byte* pComp = (byte*)(void*)bmdComp.Scan0;

            // get average pixel value of sample area
            CompAvg     = GetAverage(bmdComp, 0, bmdComp.Width, 0, bmdComp.Height); 

            for (int y = 0; y < SampleHeight; y++)
            {
                for (int x = 0; x < SampleWidth; x++)
                {
                    int Sidx = (y * bmdSample.Stride) + x * m_iPixelSize;

                    CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * (pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg);
                    a =  pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg;
                    CN2Sum += (a * a);
                }
            }

            // release pixeldata of current search area (commented out when using moving average)
            m_bmpCompositeRaw.UnlockBits(bmdComp);

            float CN2 = (float)Math.Sqrt(CN2Sum);
            float CN0 = (float)CN0Sum;
            tmpRes = CN0 / (CN1 * CN2);

            if (tmpRes > res) { res = tmpRes; SampleMatchLocation.X = m_rectSearchArea.X + j; SampleMatchLocation.Y = m_rectSearchArea.Y + i; }

            // exit early if perfect match found
            if (res == 1)
            {
                m_bmpSampleRaw.UnlockBits(bmdSample);
                CCCoefficient = res;
                return;
            }
        }
    }

    m_bmpSampleRaw.UnlockBits(bmdSample);
    CCCoefficient = res;
}
* +1036 * [/ EDIT] Метод корреляции (оригинал):
float res = 0;
float tmpRes = 0;

// get bit data of sample area
BitmapData bmdSample = m_bmpSampleRaw.LockBits(m_rectSampleArea, ImageLockMode.ReadOnly, m_bmpSampleRaw.PixelFormat);

// calculate sample average and coefficient 1 (stays same for all iterations)
int SampleAvg  = GetAverage(bmdSample, 0, bmdSample.Width, 0, bmdSample.Height);
float CN1      = GetCN1(bmdSample, SampleAvg);

int CompAvg = 0;
BitmapData bmdComp = null;
Rectangle compRect;

unsafe
{
// iterate through search area (I know it skips if areas have same size)
for (int i = 0; i < (m_rectSearchArea.Height - m_rectSampleArea.Height); i++)
{
    for (int j = 0; j < (m_rectSearchArea.Width - m_rectSampleArea.Width); j++)
    {
        int CN0Sum = 0;
        int CN2Sum = 0;

        // create composite pixel data at current search location
        compRect    = new Rectangle(m_rectSearchArea.X + j, m_rectSearchArea.Y + i,   m_rectSampleArea.Width, m_rectSampleArea.Height);
        bmdComp     = m_bmpCompositeRaw.LockBits(compRect, ImageLockMode.ReadOnly, m_bmpCompositeRaw.PixelFormat);

        CompAvg     = GetAverage(bmdComp, 0, bmdComp.Width, 0, bmdComp.Height);

        // the actual correlation loops
        byte* pSample = (byte*)(void*)bmdSample.Scan0;
        byte* pComp   = (byte*)(void*)bmdComp.Scan0;
        for (int y = 0; y < bmdSample.Height; y++)
        {
            for (int x = 0; x < bmdSample.Width; x++)
            {
                int Sidx = (y * bmdSample.Stride) + x * m_iPixelSize; // same stride assumed
                //CN0Sum += (GetPixelValue(pSample, Sidx) - SampleAvg) * (GetPixelValue(pComp, Sidx) - CompAvg);
                CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * (pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg);
                //CN2Sum += (long)Math.Pow((GetPixelValue(pComp, Sidx) - CompAvg), 2);
                CN2Sum += (int)Math.Pow((pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg), 2);
            }
        }

        // release pixeldata of current search area
        m_bmpCompositeRaw.UnlockBits(bmdComp);
        tmpRes = (float)CN0Sum / (CN1 * (float)Math.Sqrt(CN2Sum));

        if (tmpRes > res) { res = tmpRes; SampleMatchLocation.X = m_rectSearchArea.X + j; SampleMatchLocation.Y = m_rectSearchArea.Y + i; }

        // exit early if perfect match found
        if (res == 1)
        {
            m_bmpSampleRaw.UnlockBits(bmdSample);
            CCCoefficient = res;
            return;
        }
    }
}
} // unsafe
m_bmpSampleRaw.UnlockBits(bmdSample);
CCCoefficient = res;

Метод, используемый для расчета среднего значения по указанной области:

private int GetAverage(BitmapData bmpData, int X1, int X2, int Y1, int Y2)
{
    int total = 0;
    if (X2 == 0 || X2 == X1) X2++;
    if (Y2 == 0 || Y2 == Y1) Y2++;
    unsafe
    {
        byte* p = (byte*)(void*)bmpData.Scan0;
        for (int y = Y1; y < Y2; y++)
        {
            for (int x = X1; x <X2; x++)
            {
                int idx = (y * bmpData.Stride) + x * m_iPixelSize;
                //total += GetPixelValue(p, idx);
                total += p[idx] + p[idx + 1] + p[idx + 2];
            }
        }
    }
    return total / ((X2 - X1) * (Y2 - Y1));
}

Небольшая функция для вычисления средних значений по пикселям, быстро отбрасываемая:

private unsafe Int32 GetPixelValue(byte* pPixel, int idx)
{
    // add up all color values and return
    return pPixel[idx] + pPixel[idx + 1] + pPixel[idx + 2];
}

Функция, используемая для расчета неизменной части уравнения

private float GetCN1(BitmapData bmpData, long avg)
{
    double Sum = 0;
    unsafe
    {
        byte* p = (byte*)(void*)bmpData.Scan0;
        for (int y = 0; y < bmpData.Height; y++)
        {
            for (int x = 0; x < bmpData.Width; x++)
            {
                int idx = (y * bmpData.Stride) + x * m_iPixelSize;
                Sum += Math.Pow(p[idx] + p[idx + 1] + p[idx + 2] - avg, 2);
            }
        }
    }
    return (float)Math.Sqrt(Sum);
}

Ответы [ 4 ]

3 голосов
/ 23 марта 2012

О производительности и «четырех вложенных циклах»:

Вычислительная сложность корреляции, вычисленная «по определению», является произведением всех размеров изображения и рисунка O (W * H * PW * PH). Но есть быстрый метод, использующий БПФ (быстрое преобразование Фурье) со сложностью O (N ^ 2 * Log (N)), где N - наибольшее измерение.

Шаги:

Нулевое заполнение (для выравнивания размеров)

БПФ изображения и рисунка;

Комплексное размножение FT каждого компонента с комплексно-сопряженным шаблоном FT;

БПФ назад сложного продукта;

Нормализация

Дополнительно: Часто полезно сдвинуть вниз все значения в матрице - найти среднее значение и вычесть его из всего значения, чтобы получить «биполярный» сигнал. В противном случае максимальное значение корр. матрица может находиться в пиковых значениях исходной матрицы, а не в искомой позиции фрагмента

0 голосов
/ 24 марта 2012

Вы хотите избежать выполнения избыточной математики и выполнения избыточных вызовов функций, поэтому это:

for (int i = 0; i < (m_rectSearchArea.Height - m_rectSampleArea.Height); i++)

должно быть больше похоже на:

int height = m_rectSearchArea.Height - m_rectSampleArea.Height;
for (int i = 0; i < height; i++)

редактировать

Вы также можете попробовать заменить:

int Sidx = (y * bmdSample.Stride) + x * m_iPixelSize; // same stride assumed
//CN0Sum += (GetPixelValue(pSample, Sidx) - SampleAvg) * (GetPixelValue(pComp, Sidx) - CompAvg);
CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * (pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg);
//CN2Sum += (long)Math.Pow((GetPixelValue(pComp, Sidx) - CompAvg), 2);
CN2Sum += (int)Math.Pow((pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg), 2);

с:

var a = pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg; // this may already be happening
CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * a;
CN2Sum += (int)(a * a); // replacing a function call with a multiply will get you a little speed

Следует иметь в виду, что JITer не слишком хорош для такого рода кода, и вы, вероятно, сможете добиться большего, перенеся часть вашего проекта на C и P / Invoking из приложения C #.

0 голосов
/ 23 марта 2012

Как и многие алгоритмы, связанные с изображениями, похоже, что его можно легко распараллелить.

Если это правда, запуск этого на графическом процессоре даст вам огромное увеличение скорости ... если вы готовы пойти так далеко.


Вы не можете запускать код .Net непосредственно на графическом процессоре, но есть библиотеки , которые преобразуют ваш код в нечто, работающее на GPU. В противном случае вам нужно выучить язык шейдеров

0 голосов
/ 23 марта 2012

В этой области:

for (int y = 0; y < bmdSample.Height; y++)
{
   for (int x = 0; x < bmdSample.Width; x++)
   {
     int Sidx = (y * bmdSample.Stride) + x * m_iPixelSize; // same stride assumed
     CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * (pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg);
     CN2Sum += (int)Math.Pow((pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg), 2);
   }
}

Вы используете два цикла, когда с 1 все будет в порядке - поскольку Sidx колеблется от 0 до (bmdSample.Height * bmdSample.Stride) + bmdSample.Width * m_iPixelSize, у вас может быть только один цикл. Без расчета на Sidx. Это должна быть та же функциональность:

for (int Sidx = 0; Sidx < (bmdSample.Height * bmdSample.Stride) + bmdSample.Width * m_iPixelSize; Sidx++)
{
   CN0Sum += (pSample[Sidx] + pSample[Sidx + 1] + pSample[Sidx + 2] - SampleAvg) * (pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg);
   CN2Sum += (int)Math.Pow((pComp[Sidx] + pComp[Sidx + 1] + pComp[Sidx + 2] - CompAvg), 2);
}

Вы можете сделать подобный "трюк" с GetAverage и GetCN1

...