On the first article[ASIFT 1], IM_X=800 IM_Y=600;
///// Resize if the resize flag is not set or if the flag is set unequal to 0
float wS = IM_X;
float hS = IM_Y;
float zoom1=0, zoom2=0;
int wS1=0, hS1=0, wS2=0, hS2=0;
vector<float> ipixels1_zoom, ipixels2_zoom;
int flag_resize = 1;
if (argc == 9)
{
flag_resize = atoi(argv[8]);
}
Whether the images are resized or not is decided by the input parameters. If there are 9 parameters, the last one will decide the flag_resize by getting the int from sting (atoi function) .
If there are 8 parameters, resizing the image by ws and hs, namely 800*600.
if ((argc == 8) || (flag_resize != 0))
{
cout << "WARNING: The input images are resized to " << wS << "x" << hS << " for ASIFT. " << endl
<< " But the results will be normalized to the original image size." << endl << endl;
float InitSigma_aa = 1.6;
float fproj_p, fproj_bg;
char fproj_i;
float *fproj_x4, *fproj_y4;
int fproj_o;
fproj_o = 3;
fproj_p = 0;
fproj_i = 0;
fproj_bg = 0;
fproj_x4 = 0;
fproj_y4 = 0;
float areaS = wS * hS;
// Resize image 1
float area1 = w1 * h1;
zoom1 = sqrt(area1/areaS);
wS1 = (int) (w1 / zoom1);
hS1 = (int) (h1 / zoom1);
int fproj_sx = wS1;
int fproj_sy = hS1;
float fproj_x1 = 0;
float fproj_y1 = 0;
float fproj_x2 = wS1;
float fproj_y2 = 0;
float fproj_x3 = 0;
float fproj_y3 = hS1;
zoom1 shows the rate between now and past.
zoom1>1, image becomes smaller. We should change sigma_aa (smaller) and GaussianBlur1D. (反走样高斯滤波)。空间反走样滤波在图形从高分辨率变成低分辨率的过程中很常用,ASIFT使用Gauss滤波器,来做这个滤波。关于spcalial Anti-aliasing 可以查看WIKIPEDIA.
then we can get a new ipixels1.
/* Anti-aliasing filtering along vertical direction */
if ( zoom1 > 1 )
{
float sigma_aa = InitSigma_aa * zoom1 / 2;
GaussianBlur1D(ipixels1,w1,h1,sigma_aa,1);
GaussianBlur1D(ipixels1,w1,h1,sigma_aa,0);
}
if zoom1<=1. we would not change InitSigma_aa.
函数中参数的意思:ipixels1图像原始数据,作为输入; ipixels1_zoom图形resize过后的数据,作为输出
w1, h1图像原始的宽和高, &fproj_sx, &fproj_sy图像resize过后的宽和高,fproj_x2, fproj_y3图像resize过后的宽和高&fproj_bg =0, &fproj_o=3, &fproj_p=0,
&fproj_i=0 , fproj_x1=0 , fproj_y1=0 , fproj_x2 , fproj_y2=0 , fproj_x3=0 ,fproj_y3, fproj_x4=0, fproj_y4=0
// simulate a tilt: subsample the image along the vertical axis by a factor of t.
ipixels1_zoom.resize(wS1*hS1);
fproj (ipixels1, ipixels1_zoom, w1, h1, &fproj_sx, &fproj_sy, &fproj_bg, &fproj_o, &fproj_p,
&fproj_i , fproj_x1 , fproj_y1 , fproj_x2 , fproj_y2 , fproj_x3 , fproj_y3, fproj_x4, fproj_y4);
void fproj(vector<float>& in, vector<float>& out, int nx, int ny, int *sx, int *sy, float *bg, int *o, float *p, char *i, float X1, float Y1, float X2, float Y2, float X3, float Y3, float *x4, float *y4)
/* Fimage in,out;
int *sx,*sy,*o;
char *i;
float *bg,*p,X1,Y1,X2,Y2,X3,Y3,*x4,*y4; */
{
/* int n1,n2,nx,ny,x,y,xi,yi,adr,dx,dy;*/
int n1,n2,x,y,xi,yi,adr,dx,dy;
float res,xx,yy,xp,yp,ux,uy,a,b,d,fx,fy,x12,x13,y12,y13;
float cx[12],cy[12],ak[13];
/* Fimage ref,coeffs; */
// float *ref, *coeffs;
vector<float> ref, coeffs;
/* CHECK ORDER */
if (*o!=0 && *o!=1 && *o!=-3 &&
*o!=3 && *o!=5 && *o!=7 && *o!=9 && *o!=11)
/* mwerror(FATAL,1,"unrecognized interpolation order.\n"); */
{
printf("unrecognized interpolation order.\n");
exit(-1);
}
/* ALLOCATE NEW IMAGE */
/* nx = in->ncol; ny = in->nrow; */
/* out = mw_change_fimage(out,*sy,*sx);
if (!out) mwerror(FATAL,1,"not enough memory\n"); */
if (*o>=3) {
/* coeffs = mw_new_fimage();
finvspline(in,*o,coeffs); */
// coeffs = new float[nx*ny];
coeffs = vector<float>(nx*ny);
finvspline(in,*o,coeffs,nx,ny);
ref = coeffs;
if (*o>3) init_splinen(ak,*o);
} else {
// coeffs = NULL;
ref = in;
}
/* COMPUTE NEW BASIS */
if (i) {
x12 = (X2-X1)/(float)nx;
y12 = (Y2-Y1)/(float)nx;
x13 = (X3-X1)/(float)ny;
y13 = (Y3-Y1)/(float)ny;
} else {
x12 = (X2-X1)/(float)(*sx);
y12 = (Y2-Y1)/(float)(*sx);
x13 = (X3-X1)/(float)(*sy);
y13 = (Y3-Y1)/(float)(*sy);
}
if (y4) {
xx=((*x4-X1)*(Y3-Y1)-(*y4-Y1)*(X3-X1))/((X2-X1)*(Y3-Y1)-(Y2-Y1)*(X3-X1));
yy=((*x4-X1)*(Y2-Y1)-(*y4-Y1)*(X2-X1))/((X3-X1)*(Y2-Y1)-(Y3-Y1)*(X2-X1));
a = (yy-1.0)/(1.0-xx-yy);
b = (xx-1.0)/(1.0-xx-yy);
}
else
{
a=b=0.0;
}
/********** MAIN LOOP **********/
for (x=0;x<*sx;x++)
for (y=0;y<*sy;y++) {
/* COMPUTE LOCATION IN INPUT IMAGE */
if (i) {
xx = 0.5+(((float)x-X1)*y13-((float)y-Y1)*x13)/(x12*y13-y12*x13);
yy = 0.5-(((float)x-X1)*y12-((float)y-Y1)*x12)/(x12*y13-y12*x13);
d = 1.0-(a/(a+1.0))*xx/(float)nx-(b/(b+1.0))*yy/(float)ny;
xp = xx/((a+1.0)*d);
yp = yy/((b+1.0)*d);
} else {
fx = (float)x + 0.5;
fy = (float)y + 0.5;
d = a*fx/(float)(*sx)+b*fy/(float)(*sy)+1.0;
xx = (a+1.0)*fx/d;
yy = (b+1.0)*fy/d;
xp = X1 + xx*x12 + yy*x13;
yp = Y1 + xx*y12 + yy*y13;
}
/* INTERPOLATION */
if (*o==0) {
/* zero order interpolation (pixel replication) */
xi = (int)floor((double)xp);
yi = (int)floor((double)yp);
/* if (xi<0 || xi>=in->ncol || yi<0 || yi>=in->nrow)*/
if (xi<0 || xi>=nx || yi<0 || yi>=ny)
res = *bg;
else
/* res = in->gray[yi*in->ncol+xi]; */
res = in[yi*nx+xi];
} else {
/* higher order interpolations */
if (xp<0. || xp>(float)nx || yp<0. || yp>(float)ny) res=*bg;
else {
xp -= 0.5; yp -= 0.5;
xi = (int)floor((double)xp);
yi = (int)floor((double)yp);
ux = xp-(float)xi;
uy = yp-(float)yi;
switch (*o)
{
case 1: /* first order interpolation (bilinear) */
n2 = 1;
cx[0]=ux; cx[1]=1.-ux;
cy[0]=uy; cy[1]=1.-uy;
break;
case -3: /* third order interpolation (bicubic Keys' function) */
n2 = 2;
keys(cx,ux,*p);
keys(cy,uy,*p);
break;
case 3: /* spline of order 3 */
n2 = 2;
spline3(cx,ux);
spline3(cy,uy);
break;
default: /* spline of order >3 */
n2 = (1+*o)/2;
splinen(cx,ux,ak,*o);
splinen(cy,uy,ak,*o);
break;
}
res = 0.; n1 = 1-n2;
/* this test saves computation time */
if (xi+n1>=0 && xi+n2<nx && yi+n1>=0 && yi+n2<ny) {
adr = yi*nx+xi;
for (dy=n1;dy<=n2;dy++)
for (dx=n1;dx<=n2;dx++)
/* res += cy[n2-dy]*cx[n2-dx]*ref->gray[adr+nx*dy+dx];*/
res += cy[n2-dy]*cx[n2-dx]*ref[adr+nx*dy+dx];
} else
for (dy=n1;dy<=n2;dy++)
for (dx=n1;dx<=n2;dx++)
/* res += cy[n2-dy]*cx[n2-dx]*v(ref,xi+dx,yi+dy,*bg); */
res += cy[n2-dy]*cx[n2-dx]*v(ref,xi+dx,yi+dy,*bg,nx,ny);
}
}
/* out->gray[y*(*sx)+x] = res; */
out[y*(*sx)+x] = res;
}
//if (coeffs)
/* mw_delete_fimage(coeffs); */
// delete[] coeffs;
}