Menu
▾
▴
octave-dev
|
From: <Jus...@UL...> - 2004-10-18 13:55:10
|
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
|
|
From: Stefan v. d. W. <st...@su...> - 2004-10-21 09:50:10
|
Hi Justin The script works well. If you want it to be incorporated into octave-forge - add a license - fix imgrotate() octave:1> imgrotate error: `theta' undefined near line 31 column 11 error: evaluating binary operator `*' near line 31, column 17 error: evaluating binary operator `/' near line 31, column 21 error: evaluating assignment expression near line 31, column 9 error: called from `imgrotate' in file `/tmp/imgrotate.m' - attach it to an email (do not include it in the text body) Regards Stefan On Mon, Oct 18, 2004 at 03:55:03PM +0200, Justus H. Piater wrote: > 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. > > 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 > |
|
From: <Jus...@UL...> - 2004-10-26 07:51:02
Attachments:
imrotate.m
|
OK, here's my proposed replacement for imrotate, with fixes applied, humbly submitted for inclusion in octave-forge. Justus Stefan van der Walt <st...@su...> wrote on Thu, 21 Oct 2004 11:49:51 +0200: > Hi Justin > > The script works well. If you want it to be incorporated into > octave-forge |
|
From: Etienne G. <et...@cs...> - 2004-10-27 12:07:43
|
Hi all,
I have just tested imrotate briefly and commited it on octave-forge
cvs. Thanks for your contribution,
Etienne
On Tue, Oct 26, 2004 at 09:50:41AM +0200, Justus H. Piater wrote:
# OK, here's my proposed replacement for imrotate, with fixes applied,
# humbly submitted for inclusion in octave-forge.
#=20
# Justus
#=20
#=20
# Stefan van der Walt <st...@su...> wrote on Thu, 21 Oct 2004
# 11:49:51 +0200:
#=20
# > Hi Justin
# >
# > The script works well. If you want it to be incorporated into
# > octave-forge
#=20
# ## Copyright (C) 2004 Justus Piater
# ##
# ## This program is free software; you can redistribute it and/or
# ## modify it under the terms of the GNU General Public License
# ## as published by the Free Software Foundation; either version 2
# ## of the License, or (at your option) any later version.
# ##
# ## This program is distributed in the hope that it will be useful, but
# ## WITHOUT ANY WARRANTY; without even the implied warranty of
# ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# ## General Public License for more details.
# ##
# ## You should have received a copy of the GNU General Public License
# ## along with this program; if not, write to the Free Software
# ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-130=
7, USA.
#=20
# ## -*- texinfo -*-
# ## @deftypefn {Function File} {}=20
# ## imrotate(@var{imgPre}, @var{theta}, @var{method}, @var{bb=
ox})
# ## Rotation of a 2D matrix about its center.
# ##
# ## Input parameters:
# ##
# ## @var{imgPre} an input image matrix
# ##
# ## @var{theta} the rotation angle in degrees counterclockwise
# ##
# ## @var{method} "nearest" neighbor (default; faster, but produces
# ## aliasing effects) or "bilinear" interpolation
# ## (preferred; does anti-aliasing)
# ##
# ## @var{bbox} "loose" (default) or "crop"
# ##
# ## Output parameters:
# ##
# ## @var{imgPos}t the rotated image matrix
# ##
# ## @var{H} the homography mapping original to rotated pixel
# ## coordinates. To map a coordinate vector c =3D [x;y=
] to its
# ## rotated location, compute round((@var{H} * [c; 1])(1:2)).
# ## @end deftypefn
#=20
# ## Author: Justus H. Piater <Jus...@UL...>
# ## Created: 2004-10-18
# ## Version: 0.1
#=20
# function [imgPost, H] =3D imrotate(imgPre, theta, method, bbox)
# if (nargin < 2)
# help("imrotate");
# return;
# endif
#=20
# theta =3D theta * pi/180;
#=20
# sizePre =3D size(imgPre);
#=20
# ## We think in x,y coordinates here (rather than row,column).
#=20
# R =3D [cos(theta) sin(theta); -sin(theta) cos(theta)];
#=20
# if (nargin =3D=3D 4 && strcmp(bbox, "crop"))
# sizePost =3D sizePre;
# else
# ## Compute new size by projecting image corners through the rotatio=
n:
# 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
#=20
# ## 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;
#=20
# ## And here is the homography mapping old to new coordinates:
# H =3D [[R; 0 0] [T; 1]];
#=20
# Hinv =3D inv(H);
#=20
# ## "Pre" variables hold pre -rotation values;
# ## "Post" variables hold post-rotation values.
#=20
# ## Target coordinates:
# [xPost, yPost] =3D meshgrid(1:(sizePost(2)), 1:(sizePost(1)));
#=20
# ## 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].
#=20
# ## 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));
#=20
# imgPost =3D zeros(sizePost);
# imgPost(iPost) =3D imgPre(iPre);
# else
# ## bilinear interpolation between the four floor and ceiling coordi=
nates
# xPreFloor =3D floor(xPre);
# xPreCeil =3D ceil (xPre);
# yPreFloor =3D floor(yPre);
# yPreCeil =3D ceil (yPre);
#=20
# valid =3D find(1 <=3D xPreFloor & xPreCeil <=3D sizePre(2) &
# 1 <=3D yPreFloor & yPreCeil <=3D sizePre(1) );
#=20
# xPreFloor =3D xPreFloor(valid);
# xPreCeil =3D xPreCeil (valid);
# yPreFloor =3D yPreFloor(valid);
# yPreCeil =3D yPreCeil (valid);
#=20
# ## 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);
#=20
# ## We'll have to weight by the fractional part of the coordinates:
# xPreFrac =3D xPre(valid) - xPreFloor;
# yPreFrac =3D yPre(valid) - yPreFloor;
#=20
# iPost =3D sub2ind(sizePost, yPost(valid), xPost(valid));
#=20
# 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
# --=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
--=20
Etienne Grossmann ------ http://www.cs.uky.edu/~etienne
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X