Perspective Transformation

Если вбить в поисковике фразу «Перспективная трансформация» половина выдачи будет посвящено психологии. Народ интересуется личностным преображением на карантине и кормит доморощенных психологов.

Поэтому Perspective Transformation. К личностной деформации не имеет никакого отношения.

Что под эти этом понимается? Это не аффинное преобразование, это произвольная деформация прямоугольника картинки по четырем точкам с адекватным искажением содержимого.

Рис.1. Перспективное преобразование

По сути надо решить систему из 8-ми уравнений. Найти коэффициенты, сделать еще ряд телодвижений. Но здесь этого точно не будет. Во-первых, это явно просится в отдельную статью «Теории и практики», во-вторых, это медленно. Тут представлена чистая практика и готовые рецепты.

Graphics32

Признанный лидер растровых преобразований для Delphi конечно же имеет в арсенале эту трансформацию. Однако, примеры, которые во множестве представлены в составе библиотеки, не самые очевидные в плане практического применения. Поэтому, как сделать перспективную трансформацию силами Graphics32 вынес в отдельную функцию и отдельный модуль (GR32ProjectiveTransform).

Предполагаю, что все пути прописаны в Options — Library правильно. Для реализации перспективной трансформации методом Graphics32 понадобятся следующие модули: GR32, GR32_Transforms, GR32_Resamplers.

Перспективная трансформация Graphics32
function GR32PerspTransformBitmap(const ABitmap: TBitmap;
  const p1,p2,p3,p4: TPointF;
  const ABilinearFilter: Boolean = False): TBitmap;
var
  Source: TBitmap32;
  Bitmap: TBitmap32;
  Transformation: TProjectiveTransformation;
begin
  Source := nil;
  Bitmap := nil;
  Transformation := nil;

  try
    Source := TBitmap32.Create;
    Source.Assign(ABitmap);
    Bitmap := TBitmap32.Create;
    Bitmap.Assign(Source);
    Transformation := TProjectiveTransformation.Create;

    if ABilinearFilter then
      TLinearResampler.Create(Source)
    else
      TNearestResampler.Create(Source);

    Transformation.X0 := p1.X;
    Transformation.Y0 := p1.Y;
    Transformation.X1 := p2.X;
    Transformation.Y1 := p2.Y;
    Transformation.X2 := p3.X;
    Transformation.Y2 := p3.Y;
    Transformation.X3 := p4.X;
    Transformation.Y3 := p4.Y;

    Transformation.SrcRect :=
      FloatRect(0, 0, Source.Width - 1, Source.Height - 1);

    Bitmap.Clear($00000000);
    Transform(Bitmap, Source, Transformation);

    Result := TBitmap.Create;
    Result.Assign(Bitmap);

  finally
    FreeAndNil(Transformation);
    FreeAndNil(Source);
    FreeAndNil(Bitmap);
  end;
end;
[свернуть]

Параметры:

ABitmap: TBitmap - исходное изображение
p1,p2,p3,p4: TPointF - точки трансформации (см.рис.1)
ABilinearFilter: Boolean - использовать билинейный фильтр

В демонстрационном проекте представлен урезанный вариант библиотеки. Распространением библиотек этот сайт не занимается. Тут всегда представлен только необходимый для работы приложения минимум. Скачать и поставить полный комплект (исходники, компоненты, примеры) можно по этой ссылке.

Рис.2. Graphics32

OpenCV

Легендарная библиотека компьютерного зрения, обработки изображений и массы всего другого. Тащит с собой внушительный набор dll. Написана на С/С++, поэтому дополнительно включает в себя ряд обязательных библиотек от VC.

Есть библиотека для Delphi. Но для устаревшей (но не потерявшей от этого своей прелести и функциональности) версии 2.4.13. При установке следуем инструкции в той части, которая касается непосредственно Delphi.

Недостающие dll для OpenCV 2.4.13 берем тут.

Библиотеки VC можно взять из демонстрационного проекта. В проекте использованы только «нужные» библиотеки OpenCV, тащить весь набор смысла нет. Зато присутствует набор для VC. Как минимум:

  • MSVCP140D.dll
  • VCRUNTIME140D.dll
  • consrt140d.dll
  • ucrtbased.dll

Также, в проекте использованы только необходимые заголовочные файлы. Если нужна библиотека в полном составе, необходимо качать по ссылкам выше.

Снова отдельный модуль OCVProjectiveTransform и функция OCVPerspTransformBitmap с абсолютно таким же набором параметров, как и в GR32PerspTransformBitmap. Есть дополнительно еще две функции, одна для перевода стандартного для Delphi TBitmap в формат представления, понятного OpenCV PIplImage:

function BitmapToIplImage(const bitmap: TBitmap): PIplImage;

Взято практически полностью из ocv.utils. Произошло это из-за того, что следующая функция реализована в модуле плохо и пришлось ее переписывать:

function cvImage2Bitmap(img: PIplImage): TBitmap;

Как не трудно догадаться, осуществляет обратное преобразование.

Для перспективных (и не только) трансформаций необходимым и достаточным является набор из следующих модулей:

  • ocv.highgui_c
  • ocv.core_c
  • ocv.core.types_c
  • ocv.imgproc_c
  • ocv.imgproc.types_c
Перспективная трансформация OpenCV
//*******************************************************************
//  Get bitmap in perspective projection (OpenCV)
//*******************************************************************
function OCVPerspTransformBitmap(const ABitmap: TBitmap;
  const p1,p2,p3,p4: TPointF;
  const ABilinearFilter: Boolean = False): TBitmap;
var
  src: pIplImage;
  dst: pIplImage;
  srcQuad, dstQuad: pCvPoint2D32f;
  warp_matrix: pCvMat;
  Flags: Integer;
begin
  src := nil;
  dst := nil;

  try
    // prepare bitmap
    src := BitmapToIplImage(ABitmap);
    // initialize transformation matrix
    warp_matrix := cvCreateMat(3, 3, CV_32FC1);
    // clone image
    dst := cvCloneImage(src);
    // create points array
    srcQuad := AllocMem(SizeOf(TCvPoint2D32f) * 4);
    dstQuad := AllocMem(SizeOf(TCvPoint2D32f) * 4);
    // init source points
    srcQuad[0].x := 0;             // src Top left
    srcQuad[0].y := 0;
    srcQuad[1].x := src.width - 1; // src Top right
    srcQuad[1].y := 0;
    srcQuad[2].x := 0;             // src Bottom left
    srcQuad[2].y := src.height - 1;
    srcQuad[3].x := src.width - 1; // src Bottom right
    srcQuad[3].y := src.height - 1;

    // destination initialization
    dstQuad[0].x := p1.x;          // dst Top left
    dstQuad[0].y := p1.y;
    dstQuad[1].x := p2.x;          // dst Top right
    dstQuad[1].y := p2.y;
    dstQuad[2].x := p4.x;          // dst Bottom left
    dstQuad[2].y := p4.y;
    dstQuad[3].x := p3.x;          // dst Bottom right
    dstQuad[3].y := p3.y;

    // get transformation matrix
    cvGetPerspectiveTransform(srcQuad, dstQuad, warp_matrix);
    // init quality warp
    Flags := CV_WARP_FILL_OUTLIERS;
    if ABilinearFilter then
      Flags := Flags or CV_INTER_LINEAR;
    // warp
    cvWarpPerspective(src, dst, warp_matrix,
      flags, cvScalarAll(0));

    // get bitmap
    Result := cvImage2Bitmap(dst);

    // free all
    FreeMem(srcQuad);
    FreeMem(dstQuad);
    cvReleaseImage(src);
    cvReleaseImage(dst);
    cvReleaseMat(warp_matrix);

  except
    Result := TBitmap.Create;
    Result.Assign(ABitmap);
  end;
end;
[свернуть]

При написании функции OCVPerspTransformBitmap за основу был взят пример cv_WarpPerspective (samples\LibDemo\cvWrapPrespective). Безусловно, намного более очевидный, чем в Graphics32.

ВНИМАНИЕ! В демо примере используются библиотеки только для режима DEBUG (буковка d в конце названия библиотек) . Чтобы заработало в RELEASE — удалите эту буковку или продублируйте библиотеки и удалите «d». Однако, самый лучший способ, это взять полный набор библиотек, ссылка на которые представлена выше.

DEBUGRELEASE
opencv_core2413d.dllopencv_core2413.dll
opencv_highgui2413d.dllopencv_highgui2413.dll
opencv_imgproc2413d.dllopencv_imgproc2413.dll

Сравнение Graphics32 и OpenCV

Рис.3. OpenCV

Сравним рис.2. и рис.3. Два трансформированных изображения. Координаты и исходные изображения одни и те же. Разные методы получения. Сравним время:

МетодНебоГород
Graphics320.02810.0189
OpenCV0.11510.0317
Таб.1. Сравнение времени работы методов

OpenCV явно проигрывает. С чем это связано. Если заглянуть в исходники, то можно убедиться, что OpenCV «честно» решает систему линейных уравнений. По сути, решает проблему «в лоб».

Если были скачаны исходники библиотек OpenСV 2.4.13, то в модуле opencv\sources\modules\imgproc\src\imgwarp.cpp находим такой код:

Составление матрицы уравнений перспективной трансформации в OpenCV
/* Calculates coefficients of perspective transformation
 * which maps (xi,yi) to (ui,vi), (i=1,2,3,4):
 *
 *      c00*xi + c01*yi + c02
 * ui = ---------------------
 *      c20*xi + c21*yi + c22
 *
 *      c10*xi + c11*yi + c12
 * vi = ---------------------
 *      c20*xi + c21*yi + c22
 *
 * Coefficients are calculated by solving linear system:
 * / x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
 * | x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
 * | x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
 * | x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,
 * |  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
 * |  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
 * |  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
 * \  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/
 *
 * where:
 *   cij - matrix coefficients, c22 = 1
 */
cv::Mat cv::getPerspectiveTransform( const Point2f src[], 
  const Point2f dst[] )
{
    Mat M(3, 3, CV_64F), X(8, 1, CV_64F, M.data);
    double a[8][8], b[8];
    Mat A(8, 8, CV_64F, a), B(8, 1, CV_64F, b);

    for( int i = 0; i < 4; ++i )
    {
        a[i][0] = a[i+4][3] = src[i].x;
        a[i][1] = a[i+4][4] = src[i].y;
        a[i][2] = a[i+4][5] = 1;
        a[i][3] = a[i][4] = a[i][5] =
        a[i+4][0] = a[i+4][1] = a[i+4][2] = 0;
        a[i][6] = -src[i].x*dst[i].x;
        a[i][7] = -src[i].y*dst[i].x;
        a[i+4][6] = -src[i].x*dst[i].y;
        a[i+4][7] = -src[i].y*dst[i].y;
        b[i] = dst[i].x;
        b[i+4] = dst[i].y;
    }

    solve( A, B, X, DECOMP_SVD );
    ((double*)M.data)[8] = 1.;

    return M;
}
[свернуть]

Что имеем: с одной стороны очень быстрый Graphics32, которые требует распространяться под такой же лицензией и открытым кодом, с другой — более медленный OpenCV, который всего этого не требует, но тащит за собой внушительный прицеп из dll.

Однако, если использовать модули из указанной выше библиотеки OpenCV для Delphi, то помимо прицепа c dll, действуют те же правила, что и для Graphics32 — и код открытый, и лицензия такая же.

Конечно пишем свой вариант. Зачем? Есть ряд вещей (о которых частично ниже), которые меня не устраивают в реализации. Второе, не устраивает лицензия.

IP76

Составление матрицы системы уравнений, ее решение, нахождение коэффициентов — путь явно более медленный, чем тот, который предлагает Graphics32 (см.таб.1).

Поэтому за основу взял алгоритм нахождения матрицы трансформации Graphics32. Алгоритм сильно упрощен, переделан, лишен всех наворотов с целочисленной арифметикой. Без применения билинейного фильтра работает быстрее, чем алгоритм Graphics32 (конечно, также без билинейной растеризации).

Рис.4. IP76 vs Graphics32

На рисунке плохо видно, но IP76: 10.6 миллисекунд, GR32: 11.8 миллисекунд. Чем больше картинка, тем выигрыш ощутимей.

Рис.5. IP76 vs Graphics32 (7.9 msec / 11.4 msec)
Матрица перспективной трансформаций IP76
//*******************************************************************
//  get transorm matrix
//*******************************************************************
function GetPerspTransformMatrix(const p1,p2,p3,p4: TPointF;
  const ARect: TRectF): TMatrix3F;
var
  dx1, dx2, px: single;
  dy1, dy2, py: single;
  g, h, k: single;
  IM: TMatrix3F;
begin
  FillChar(Result, SizeOf(Result), 0);
  if ((abs(ARect.Width) < PrecZero) or
      (abs(ARect.Height) < PrecZero))
  then
    Exit;

  px  := p1.x - p2.x + p3.x - p4.x;
  py  := p1.y - p2.y + p3.y - p4.y;
  if (px = 0) and (py = 0) then
  begin
    Result[0, 0] := p2.x - p1.x;
    Result[1, 0] := p3.x - p2.x;
    Result[2, 0] := p1.x;

    Result[0, 1] := p2.y - p1.y;
    Result[1, 1] := p3.y - p2.y;
    Result[2, 1] := p1.y;

    Result[0, 2] := 0;
    Result[1, 2] := 0;
    Result[2, 2] := 1;
  end
  else
  begin
    dx1 := p2.x - p3.x;
    dx2 := p4.x - p3.x;
    dy1 := p2.y - p3.y;
    dy2 := p4.y - p3.y;

    k := dx1 * dy2 - dx2 * dy1;
    if abs(k) > PrecZero then
    begin
      g := (px * dy2 - py * dx2) / k;
      h := (py * dx1 - px * dy1) / k;

      Result[0, 0] := p2.x - p1.x + g * p2.x;
      Result[1, 0] := p4.x - p1.x + h * p4.x;
      Result[2, 0] := p1.x;

      Result[0, 1] := p2.y - p1.y + g * p2.y;
      Result[1, 1] := p4.y - p1.y + h * p4.y;
      Result[2, 1] := p1.y;

      Result[0, 2] := g;
      Result[1, 2] := h;
      Result[2, 2] := 1;
    end
  end;

  IM := IdentityMatrix;
  IM[0, 0] := 1 / ARect.Width;
  IM[1, 1] := 1 / ARect.Height;
  Result := MultMx(Result, IM);

  IM := IdentityMatrix;
  IM[2, 0] := -ARect.Left;
  IM[2, 1] := -ARect.Top;
  Result := MultMx(Result, IM);
end;
[свернуть]

Билинейный фильтр в моей реализации сделан исключительно по теории, без целочисленной оптимизации, что чуть влияет на скорость.

Билинейный фильтр
//**************************************************************
//  calculate color for 4 pixels (bilinear filter)
//**************************************************************
function GetBilinearColor(const ASrc: PByte; const u, v: Single;
  const w, h: Integer): Integer;
const
  C = 4;

  function GetRGB(AX,AY: Integer): PRGBQuad;
  begin
    if AX < 0 then AX := 0;
    if AX >= W then AX := W-1;
    if AY < 0 then AY := 0;
    if AY >= H then AY := H-1;
    Result := PRGBQuad(Integer(ASrc) - AY * W * C + AX * C);
  end;

  function S2B (const S: Single): Byte;
  begin
    if s < 0 then Result := 0
    else if s > 254 then Result := 255
    else Result := Round(S);
  end;

var
  x,y: Integer;
  ur, vr: Single;
  ul, vl: Single;
  v1, v2, v3, v4: Single;
  r,g,b,a: Single;
  p1, p2, p3, p4: PRGBQuad;
begin
  x := Trunc(u);
  y := Trunc(v);

  ur := u - x;
  vr := v - y;
  ul := 1 - ur;
  vl := 1 - vr;

  v1 := ul*vl;
  v2 := ur*vl;
  v3 := ul*vr;
  v4 := ur*vr;

  p1 := GetRGB(x,y);
  p2 := GetRGB(x+1,y);
  p3 := GetRGB(x,y+1);
  p4 := GetRGB(x+1,y+1);

  r := p1^.rgbRed * v1  + p2^.rgbRed * v2 +
       p3^.rgbRed * v3  + p4^.rgbRed * v4;
  g := p1^.rgbGreen * v1  + p2^.rgbGreen * v2 +
       p3^.rgbGreen * v3  + p4^.rgbGreen * v4;
  b := p1^.rgbBlue * v1  + p2^.rgbBlue * v2 +
       p3^.rgbBlue * v3  + p4^.rgbBlue * v4;

  // that is more correct:
    a := p1^.rgbReserved * v1  + p2^.rgbReserved * v2 +
         p3^.rgbReserved * v3  + p4^.rgbReserved * v4;

  // it's faster, but there may be border issues:
  // a := p1^.rgbReserved;

  Result := S2B(b) or (S2B(g) shl 8) or (S2B(r) shl 16) or 
   (S2B(a) shl 24);
end;

[свернуть]

Функция перспективной трансформации:

Создание изображения в перспективной трансформации IP76
//*******************************************************************
//  create bitmap in perspective transformation
//*******************************************************************
function GetPerspTransformBitmap(const ABitmap: TBitmap;
  const p1,p2,p3,p4: TPointF;
  const ABilinearFilter: Boolean = False): TBitmap;
const
  C = 4;
var
  GetIntColor: TGetIntColor;
  RMatrix: TMatrix3F;
  IMatrix: TMatrix3F;
  s, d, p: PInteger;
  dw, dh: Integer;
  w, h: Integer;
  x, y: Integer;
  dst: TRect;
  pnt: TPointF;
  clr: Integer;
begin
  ABitmap.PixelFormat := pf32bit;
  w := ABitmap.Width;
  h := ABitmap.Height;

  dst := ToRect(CalcBounds(p1,p2,p3,p4));
  dw := dst.Left + dst.Width;
  dh := dst.Top + dst.Height;
  if dw < w then dw := w;
  if dh < h then dh := h;

  RMatrix := GetPerspTransformMatrix(p1, p2, p3, p4,
    RectF(0, 0, w-1, h-1));
  IMatrix := InvertMx(RMatrix);

  Result := CreateBitmap(dw, dh, pf32bit);
  dst := IntersectRect(Rect(0, 0, dw, dh), dst);

  s := ABitmap.ScanLine[0];
  d := Result.ScanLine[0];

  if (dst.Width<=0) or (dst.Height<=0) then
  begin
    Result.Assign(ABitmap);
    Exit;
  end;

  if ABilinearFilter then
    GetIntColor := GetBilinearColor
  else
    GetIntColor := GetNearestColor;

  for y := dst.Top to dst.Bottom - 1 do
  begin
    p := PInteger(Integer(d) - Y * dw * C + dst.Left * C);

    for x := dst.Left to dst.Right - 1 do
    begin
      pnt := CalcSrcCoord(x, y, IMatrix);
      if (pnt.X >= 0) and (pnt.X < w) and
         (pnt.Y >= 0) and (pnt.Y < h)
      then
        clr := GetIntColor(PByte(s), pnt.x, pnt.y, w, h)
      else
        clr := 0;

      p^ := clr;
      Inc(p);
    end;
  end;
end;
[свернуть]

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

Рис.6. IP76

Расширим таблицу 1 с учетом своей реализации:

МетодНебоГород
Graphics320.02810.0189
OpenCV0.1151
(4,1 * GR32)
0.0317
(1.68 * GR32)
IP760.0364
(1.3 * GR32)
0.0273
(1.44 * GR32)
Таб.2. Сравнение времени работы всех методов

Своя реализация с применением билинейной фильтра примерно в 1.4 медленнее чем Graphics32 и сильно быстрее, чем OpenCV.

Глюки реализаций

Алгоритм и Graphics32, и OpenCV, построен таким образом, что если фигура, отображаемая на экране имеет размеры большие, чем исходное изображение, вывод на экран «обрезается».

Рис.7. Graphics32 (22.8 msec)
Рис.8. OpenCV (270.5 msec)

Проблему можно решить увеличив исходное изображение. Но это как-то искусственно и «костыльно». Либо поправив руками исходники, что тоже не самый лучший путь.

Конечно, эта проблема решена в IP76.

Рис.9. IP76 (34.7 msec)

«Шедевры»

Рис.10. Река-небо, небо-отражение заката в реке
Рис.11. Потерянный

Коллеги, спасибо за внимание!

Информация о новых статьях всегда и только в телеграм-канале. Никаких рассылок на почту.

Не забываем комментировать, ругаться и подписываться )))


Скачать (4.3 Мб): Исходники (Delphi XE 7-10)

Скачать (7.8 Мб): Исполняемый файл


0 0 голос
Рейтинг статьи
Подписаться
Уведомить о
guest
0 комментариев
Межтекстовые Отзывы
Посмотреть все комментарии
0
Оставьте комментарий! Напишите, что думаете по поводу статьи.x
()
x