桶形畸变的校正方法

%%% ================================================================%%%
%%% This arithmetic is programmed by Xiang Li
%%% Copyright 2005All Copyright Reserved by Xiang Li
%%% You can modify and distibute this programme under the GPL license
%%% If you have any question about this programme
%%% please contact with XiangLi@uestc.edu.cn
%%% ================================================================%%%


%%% ================================================================%%%
%%% first step:
%%% assume we have the inverse alias functionthen we need to determine and
%%% establish the range of x and y of the original image from this function
%%% ================================================================%%%

%%% read alias image into current working memory space
origImg = imread（‘lena_barrel.tif‘‘tif‘）;
figure;
imshow（origImg‘truesize‘）;
title（‘barrel image‘）;

%%% get width and height of alias image
origWidth = size（origImg1）;
origHeight = size（origImg2）;

%%% calculate the middle coordinate of alias image
origMidOfX = round（origWidth/2）;
origMidOfY = round（origHeight/2）;

%%% shift the coordinate to the center of the alias image
[xOrigIndexyOrigIndex] = meshgrid（1:origWidth1:origHeight）;
xNewIndex = xOrigIndex（:） - origMidOfX;
yNewIndex = yOrigIndex（:） - origMidOfY;

%%% transform rectangular to polar coordinatenote that theta is measured
%%% in radiansand please note that we can think Cartesian ordinate is just
%%% the rectangular coordinatewe can of course think it like thisanyway
[thetaradius] = cart2pol（xNewIndexyNewIndex）;

%%% a1 and a2 are image alias coeffientsso here we produce the aliasing image‘s
%%% radius in polar coordinates.
%%% noteinverse formula:newRadius = （a1 - sqrt（a1 ^ 2 - 4 * a2 * radius））/（2 * a2）
%%% is from formular:
%%% aliasRadius = a1 * radius - a2 * radius .^ 2
%%% which is the alias function I guessedand it seems true.
a1 = 2.0;
a2 = 0.005;
newRadius = （a1 - sqrt（a1 ^ 2 - 4 * a2 * radius））/（2 * a2）;
newTheta = theta;

%%% now we can map the aliasing rectangular coordinates to its polar
%%% counterpart
[xAliasIndexyAliasIndex] = pol2cart（newThetanewRadius）;

%%% here calculate the min and max position
xMinAliasIndex = double（round（min（xAliasIndex）））;
xMaxAliasIndex = double（round（max（xAliasIndex）））;
yMinAliasIndex = double（round（min（yAliasIndex）））;
yMaxAliasIndex = double（round（max（yAliasIndex）））;

%%% calculate the to-be-fixed image‘s width and height according to the min and max position
newWidth = xMaxAliasIndex - xMinAliasIndex + 1;
newHeight = yMaxAliasIndex - yMinAliasIndex + 1;

%%% ================================================================%%%
%%% end first step
%%% ================================================================%%%


%%% ================================================================%%%
%%% second step:
%%% for each point in the to-be-fixed imagecompute the corresponding position
%%% in th