[OctDev] New imrotate

[<Jus...@UL...>]

Hi,

I've had some problems with octave-forge's imrotate.m (lack of
precision, global gray-level alteration, inability to recover the
rotation mapping), so I wrote my own from scratch.=20

Highlights:

- 100% precise (unless you find bugs...), as it works by applying a
  rotation homography to each pixel coordinate. The homography is
  returned to the caller if desired.

- "nearest"-neighbor mapping or "bilinear" interpolation

- "loose" or "crop" bounding box

- faster than the existing imrotate (on the few informal tests I ran)

- does not rely on external functions not in core octave

- code is relatively compact and intuitive (to me at least...)

Caveats:

- Only a single-band image matrix is accepted as input for now.

- Its memory usage can still be improved. It makes heavy use of
  exhaustive index lists (to take advantage of vector notation); the
  current implementation keeps more of them around simultaneously than
  necessary.

The function is currently called imgrotate to allow it to coexist with
the existing imrotate. However, unless bugs show up, I don't see why
it should not eventually take its place.

You feedback is more than welcome.

Enjoy,
Justus


function [imgPost, H] =3D imgrotate(imgPre, theta, method, bbox)
  ## imgrotate: precise image rotation by homography mapping
  ##
  ## Usage:
  ##
  ## [imgPost, H] =3D imgrotate(imgPre, theta, method, bbox)
  ##
  ## Input parameters:
  ##
  ##   imgPre   an input image matrix
  ##
  ##   theta    the rotation angle in degrees counterclockwise
  ##
  ##   method   "nearest" neighbor (default; faster, but produces
  ##            aliasing effects) or "bilinear" interpolation
  ##            (preferred; does anti-aliasing)
  ##
  ##   bbox     "loose" (default) or "crop"
  ##
  ## Output parameters:
  ##
  ##   imgPost  the rotated image matrix
  ##
  ##   H        the homography mapping original to rotated pixel
  ##            coordinates. To map a coordinate vector c =3D [x;y] to its
  ##            rotated location, compute round((H * [c; 1])(1:2)).
=20=20
  ## Author: Justus H. Piater  <Jus...@UL...>
  ## Created: 2004-10-18
  ## Version: 0.1

  theta =3D theta * pi/180;

  sizePre  =3D size(imgPre);

  ## We think in x,y coordinates here (rather than row,column).

  R =3D [cos(theta) sin(theta); -sin(theta) cos(theta)];

  if (nargin =3D=3D 4 && strcmp(bbox, "crop"))
    sizePost =3D sizePre;
  else
    ## Compute new size by projecting image corners through the rotation:
    corners =3D [0, 0;
	       (R * [sizePre(2); 1])';
	       (R * [sizePre(2); sizePre(1)])';
	       (R * [1; sizePre(1)])'          ];
    sizePost(2) =3D ceil(max(corners(:,1))) - floor(min(corners(:,1)));
    sizePost(1) =3D ceil(max(corners(:,2))) - floor(min(corners(:,2)));
  endif

  ## Compute the translation part of the homography:
  oPre  =3D ([sizePre(2) ; sizePre(1) ] + 1) / 2;
  oPost =3D ([sizePost(2); sizePost(1)] + 1) / 2;
  T =3D oPost - R * oPre;

  ## And here is the homography mapping old to new coordinates:
  H =3D [[R; 0 0] [T; 1]];

  Hinv =3D inv(H);

  ## "Pre"  variables hold pre -rotation values;
  ## "Post" variables hold post-rotation values.

  ## Target coordinates:
  [xPost, yPost] =3D meshgrid(1:(sizePost(2)), 1:(sizePost(1)));

  ## Compute corresponding source coordinates:
  xPre =3D Hinv(1,1) * xPost + Hinv(1,2) * yPost + Hinv(1,3);
  yPre =3D Hinv(2,1) * xPost + Hinv(2,2) * yPost + Hinv(2,3);
  ## zPre is guaranteed to be 1, since the last row of H (and thus of
  ## Hinv) is [0 0 1].

  ## Now map the image, either by nearest neighbor or by bilinear
  ## interpolation:
  if (nargin < 3 || !size(method) || strcmp(method, "nearest"))
    ## nearest-neighbor: simply round Pre coordinates
    xPre =3D round(xPre);
    yPre =3D round(yPre);
    valid =3D find(1 <=3D xPre & xPre <=3D sizePre(2) &
		 1 <=3D yPre & yPre <=3D sizePre(1)  );
    iPre  =3D sub2ind(sizePre , yPre (valid), xPre (valid));
    iPost =3D sub2ind(sizePost, yPost(valid), xPost(valid));

    imgPost =3D zeros(sizePost);
    imgPost(iPost) =3D imgPre(iPre);
  else
    ## bilinear interpolation between the four floor and ceiling coordinates
    xPreFloor =3D floor(xPre);
    xPreCeil  =3D ceil (xPre);
    yPreFloor =3D floor(yPre);
    yPreCeil  =3D ceil (yPre);

    valid =3D find(1 <=3D xPreFloor & xPreCeil <=3D sizePre(2) &
		 1 <=3D yPreFloor & yPreCeil <=3D sizePre(1)  );

    xPreFloor =3D xPreFloor(valid);
    xPreCeil  =3D xPreCeil (valid);
    yPreFloor =3D yPreFloor(valid);
    yPreCeil  =3D yPreCeil (valid);

    ## In the following, FC =3D floor(x), ceil(y), etc.
    iPreFF =3D sub2ind(sizePre, yPreFloor, xPreFloor);
    iPreCF =3D sub2ind(sizePre, yPreFloor, xPreCeil );
    iPreCC =3D sub2ind(sizePre, yPreCeil , xPreCeil );
    iPreFC =3D sub2ind(sizePre, yPreCeil , xPreFloor);

    ## We'll have to weight by the fractional part of the coordinates:
    xPreFrac =3D xPre(valid) - xPreFloor;
    yPreFrac =3D yPre(valid) - yPreFloor;

    iPost =3D sub2ind(sizePost, yPost(valid), xPost(valid));

    imgPost =3D zeros(sizePost);
    imgPost(iPost) =3D ...
	round(imgPre(iPreFF) .* (1 - xPreFrac) .* (1 - yPreFrac) + ...
	      imgPre(iPreCF) .*      xPreFrac  .* (1 - yPreFrac) + ...
	      imgPre(iPreCC) .*      xPreFrac  .*      yPreFrac  + ...
	      imgPre(iPreFC) .* (1 - xPreFrac) .*      yPreFrac       );
  endif
endfunction


--=20
Justus H. Piater, Ph.D.         http://www.montefiore.ulg.ac.be/~piater/
Institut Montefiore, B28        Phone: +32-4-366-2279
Universit=E9 de Li=E8ge, Belgium    Fax:   +32-4-366-2620