[Из песочницы] OpenCV и иллюзия кругов на воде

Предлагаю читателям «Хабрахабра» статью о том, как школьная физика и OpenCV позволяют сделать иллюзию волн на воде. Основная сложность состоит в выборе красивого уравнения и переносе эффекта преломления света на границе раздела сред из трехмерного мира на плоскую картинку.
Решение задачи можно разделить на два этапа:

  • Создать карту распространения кругов/волн по воде;
  • Совместить полученную карту с заданным изображением.

Карта волн


Полагаю, из школьного курса физики все помнят, что круги на воде есть ни что иное, как распространение колебаний по поверхности воды. Строгий подход к решению задачи подразумевает решение дифференциального волнового уравнения в двухмерном пространстве, однако полная физическая достоверность в нашем случае не требуется — можно взять готовое решение из школьного учебника по физике или придумать свое уравнение.

04dbd627befe4a2db89f91f5c43b45a1.gif

где c — коэффициент затухания; L — длина волны; x0, y0 — координаты центра волны; t- время. Величину длины волны лучше изменять в зависимости от размеров картинки.

Из школьного курса помним, что косинус принимает значения [-1;1], однако интенсивность цвета [0;255], следовательно необходимо немного модифицировать нашу волновую функцию. Окончательно, яркость каждого пикселя на слое с волной задается следующим образом:

65600e2d0b04469ead6dc92b34b8045c.gif

poolDepth — толщина слоя воды, в данном случае, еще и амплитуда волны в начальный момент времени. Слишком большое значение приведет к существенным искажениям изображения и будет некрасиво, слишком маленькое — почти не повлияет на картинку. Для картинки 100×100 подошло значение poolDepth=20, для картинки 1600×1600 оказалось возможным выставить предельное значение poolDepth=127.

Волновая карта без lookup таблицы
void makeWaveMap(Mat& image)
{
   float simulPeriod = 10.0;     // Period of simulation
   static float time = 0.0;
   const float dt = 0.05;        // Time step
   float poolDepth = 20.0;
   int maxImageSize = image.cols > image.rows ? image.cols : image.rows;

   for (int i = 0; i < image.rows; i++) {
      for (int j = 0; j < image.cols; j++) {
         float radius = sqrt((i - image.rows/2)*(i - image.rows/2) + \
                             (j - image.cols/2)*(j - image.cols/2));
         float z = (1.0 + waveFunction(radius, time, maxImageSize))*poolDepth;
         image.at(i, j) = saturate_cast(z);
      }
   }

   time+= dt;
   time*= (time < simulPeriod);
}



Заметим, что размер волновой карты соответствует размеру исходного изображения. С учетом того, что волновая функция вычисляется в каждой точке получим гигантский объем вычислений уже для сравнительно небольшого изображения (в каждой точке вычисляем один квадратный корень, одну экспоненту и два косинуса — редкий ноутбук потянет картинку в 2МП при таком раскладе). На самом деле, в каждой точке явно вычислять и не нужно — уравнение ведь одномерное! На каждом временном шаге создаем таблицу значений волновой функции, а далее — линейная интерполяция.

Волновая карта с lookup таблицей
void makeWaveMapLUT(Mat& image)
{
   float simulPeriod = 10.0;     // Period of simulation
   static float time = 0.0;
   const float dt = 0.05;        // Time step
   float poolDepth = 20.0;
   int nLUT = image.cols > image.rows ? image.cols : image.rows;
   int maxImageSize = nLUT;
   float waveFuncLUT[nLUT];

   for (int i = 0; i < nLUT; i++) {
      float radius = saturate_cast(i);
      waveFuncLUT[i] = waveFunction(radius, time, maxImageSize);
   }

   for (int i = 0; i < image.rows; i++) {
      for (int j = 0; j < image.cols; j++) {
         float radius = sqrt((i - image.rows/2)*(i - image.rows/2) + \
                             (j - image.cols/2)*(j - image.cols/2));
         int iRad = cvRound(radius);
         float dR = radius - saturate_cast(iRad);
         float wF = waveFuncLUT[iRad] + (waveFuncLUT[iRad+1] - waveFuncLUT[iRad])*dR;
         float z = (1.0 + wF)*poolDepth;
         image.at(i, j) = saturate_cast(z);
      }
   }

   time+= dt;
   time*= (time < simulPeriod);
}



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

Наложение волновой карты на исходное изображение


Искажение изображения находящегося под водой происходит из-за преломления света на границе раздела. Фактически требуется решить школьную задачку о преломлении света на клине.

При нормальном падении света на ровную горизонтальную границу раздела сред искажения не происходит (изменения цвета картинки в результате интерференции здесь не рассматриваем — считаем, что лужа не глубокая да и на 100% физическую достоверность не претендуем).

520653fbe44947aa8b2eacfb6aec106c.gif

На рисунке изображено что происходит с лучом света падающим на наклонную поверхность, необходимо найти насколько сместился выходящий луч относительно входящего. Расстояние CD и есть искомое смещение, определяемое как (решение здесь):

3f73bad66fe3437ca8b5b18670d99451.gif

где n1=1 и n2=1.33 — показатели преломления первой (воздуха) и второй (воды) среды, соответственно.

Применительно к нашим волнам и пикселям уравнение можно записать как:

2a0b96a91b1c429a9e2f18feb1f6ac90.gif

В данном случае, alpha есть угол наклона касательной к волне в точке [i, j], который вычисляется как:

536a71f68b73407092ba6e0d52e17aac.gif

Разумеется, производную можно вычислить и аналитически, что повысит точность, однако это сделает трудоемким процесс изменения уравнения волны.

В данном случае, опять можем либо вычислять смещение для каждой точки изображения отдельно (void makeWaveMap (Mat&)), но лучше сделать lookup таблицу при первом вызове подпрограммы и использовать ее в дальнейшем (void makeWaveMapLUT (Mat&)).

Наложение без lookup таблицы
void blendWaveAndImage(Mat& sourceImage, Mat& targetImage, Mat& waveMap)
{
   static float rFactor = 1.33; // refraction factor of water

   for (int i = 1; i < sourceImage.rows-1; i++) {
      for (int j = 1; j < sourceImage.cols-1; j++) {
         float alpha, beta;

         float xDiff = waveMap.at(i+1, j) - waveMap.at(i, j);
         float yDiff = waveMap.at(i, j+1) - waveMap.at(i, j);

         alpha = atan(xDiff);
         beta = asin(sin(alpha)/rFactor);
         int xDisplace = cvRound(tan(alpha - beta)*waveMap.at(i, j));

         alpha = atan(yDiff);
         beta = asin(sin(alpha)/rFactor);
         int yDisplace = cvRound(tan(alpha - beta)*waveMap.at(i, j));

         Vec3b Intensity = sourceImage.at(i,j);

         /* Check whether displacement fits the image size */
         int dispNi = i + xDisplace;
         int dispNj = j + yDisplace;
         dispNi = (dispNi > sourceImage.rows || dispNi < 0 ? i : dispNi);
         dispNj = (dispNj > sourceImage.cols || dispNj < 0 ? j : dispNj);

         Intensity = sourceImage.at(dispNi, dispNj);

         targetImage.at(i,j) = Intensity;
      }
   }
}



Наложение с lookup таблицей
void blendWaveAndImageLUT(Mat& sourceImage, Mat& targetImage, Mat& waveMap)
{
   static float rFactor = 1.33; // refraction factor of water
   static float dispLUT[512];      //Lookup table for displacement
   static int nDispPoint = 512;

   for (int i = 0; i < nDispPoint; i++) {
      float diff = saturate_cast(i - 255);
      float alpha = atan(diff);
      float beta = asin(sin(alpha)/rFactor);
      dispLUT[i] =  tan(alpha - beta);
   }
   nDispPoint = 0;

   for (int i = 1; i < sourceImage.rows-1; i++) {
      for (int j = 1; j < sourceImage.cols-1; j++) {
         int xDiff = waveMap.at(i+1, j) - waveMap.at(i, j);
         int yDiff = waveMap.at(i, j+1) - waveMap.at(i, j);

         int xDisplace = cvRound(dispLUT[xDiff+255]*waveMap.at(i, j));

         int yDisplace = cvRound(dispLUT[yDiff+255]*waveMap.at(i, j));

         Vec3b Intensity = sourceImage.at(i,j);

         /* Check whether displacement fits the image size */
         int dispNi = i + xDisplace;
         int dispNj = j + yDisplace;
         dispNi = (dispNi > sourceImage.rows || dispNi < 0 ? i : dispNi);
         dispNj = (dispNj > sourceImage.cols || dispNj < 0 ? j : dispNj);

         Intensity = sourceImage.at(dispNi, dispNj);

         targetImage.at(i,j) = Intensity;
      }
   }
}



Полный текст программы, а также тестовые изображения здесь.

Результаты работы программы в картинках и цифрах


Среднее время обработки (создание карты волны и наложение на исходное изображение) одного кадра для картинки размером 1600×1600 составило 0.137cек и 0.415сек c lookup таблицами и без них, соответственно (на видео картинка 100×100 пикселей).

© Habrahabr.ru