/* Author : Paul Ortyl
* License: GPL
* Version: 0.02
* Date : 01.08.2005
* Description: estimates rotation angle of scanned image
* supports only PGM raw images (P5)
* reads stdin writes to stdout
* enables "automatic page rotation" of scanned text documents
* Requires: FFTW library (http://www.fftw.org)
* (-lfftw3 -lm)
* May be compiled with:
* cc -O3 -DNDEBUG -Wall -s -lfftw3 -lm pgmangle.c -o pgmangle
*
* Note : the algorithm is buggy, but for big images gives relatively
* good (or even perfect) estimation
* if you think, that smaller fragment of picture is enough -- use
* "pnmcut x y width height" command before passing it to pgmangle
* Bug/wish reports : mailto:ortylp@3miasto.net.pl
*/
#include
#include
#include
#include
#include
#include
#include
#define BUFSIZE (1<<12)
#define PI (3.14159265358979311)
static double angle(double *pict, int X, int Y, int Xpix, int Ypix,
double astart, double aend, int anum);
static double getangle(unsigned char *pix, int xpix, int ypix,
int xsize, int ysize,
int x0cut, int y0cut,
int x, int y, int freepix);
static void win_init(unsigned int Nx, unsigned int Ny);
static inline double win(unsigned int x,unsigned int y);
static void win_destroy(void);
static inline double asinc(double d, int N);
static void usage(void);
static void usage() {
puts("\n"
"(C) 2001-2005 Paul Ortyl \n"
"Licence: GPL\n\n"
"Usage: pgmangle [-h|--help]\n"
" -h --help\tprint this message\n\n"
"Description: pgmangle reads a PGM (P5) bitmap from standard input,\n"
"\testimates the rotation angle of the bitmap according to the pattern\n"
"\tfound in the bitmap and writes the value in degrees of the needed correction,\n"
"\twhich may be fed as a parameter to pnmrotate.\n"
"\tThe primary use of this tool is to estimate rotation angle of scanned\n"
"\ttext pages and correction of this rotation using pnmrotate.\n\n"
"Usage examples:\n"
" correct rotation of a grayscale bitmap\n"
" # pnmrotate `pgmangle corrected_scan.pnm\n"
" correct rotation of a colour bitmap\n"
" # pnmrotate `ppmtopgm corrected_scan.pnm "
);
return;
}
int main(int argc, char **argv)
{
char buf[BUFSIZE];
unsigned char * pix;
unsigned int x=0,y=0,c=0,xpix,ypix,
x0cut=0,xsize=0,y0cut=0,ysize=0,x2;
double ang;
if (argc>1) {
usage();
exit(1);
}
fgets(buf,BUFSIZE,stdin);
if ((buf[0]!='P') || (buf[1]!='5') || (buf[2]!='\n'))
{fprintf(stderr,"err\n");return 1;}
do {
fgets(buf,BUFSIZE,stdin);
} while (*buf=='#');
if (sscanf(buf,"%u%u",&xpix,&ypix)!=2)
{fprintf(stderr,"err, scanf1\n");return 1;}
fgets(buf,BUFSIZE,stdin);
if (sscanf(buf,"%u",&c)!=1)
{fprintf(stderr,"err, scanf2\n");return 1;}
if (c>255)
{fprintf(stderr,"err, too many colours\n");return 1;}
if ((pix=malloc(xpix*ypix))==NULL)
{fprintf(stderr,"err, malloc\n");return 1;}
fprintf(stderr,"picture size: x=%d\ty=%d\n",xpix,ypix);
xsize=xpix;
ysize=ypix;
x= 1<< ((unsigned int) (floor(log(xsize-1)/log(2))+1));
if (5*x/8>xsize) {x=5*x/8;}
else
if (3*x/4>=xsize) {x=3*x/4;}
else
if (7*x/8>=xsize) {x=7*x/8;}
y= 1<< ((unsigned int) (floor(log(ysize-1)/log(2))+1));
if (5*y/8>=ysize) {y=5*y/8;}
else
if (3*y/4>=ysize) {y=3*y/4;}
else
if (7*y/8>=ysize) {y=7*y/8;}
x2=2*(x/2+1);
fprintf(stderr,"zeropadding: x=%d\ty=%d\n",x,y);
#ifndef NDEBUG
fprintf(stderr,"memory: %.2f + %.2f = %.2fMB\n",xpix*ypix/1024.0/1024.0,
y*x2*sizeof(double)/1024.0/1024.0,
(xpix*ypix+y*x2*sizeof(double))/1024.0/1024.0);
#endif
{
struct sysinfo info;
unsigned long Memlimit;
sysinfo(&info);
Memlimit=info.totalram*info.mem_unit/4*3;
if (xpix*ypix+y*x2*sizeof(double)>Memlimit)
{
unsigned long int memlimit=(Memlimit-xpix*ypix)/sizeof(double);
unsigned long int square=floor(sqrt(memlimit));
unsigned long int square2;
fprintf(stderr,"too big, memlimit=0.75*RAM=%.2fMB\n",
Memlimit/1024.0/1024.0);
if ((xpix>square)&&(ypix>square))
{
square2= 1<< ((unsigned int) (floor(log(square-1)/log(2))+1));
if (7*square2/8xpix) {x=5*x/8;}
else
if (3*x/4>=xpix) {x=3*x/4;}
else
if (7*x/8>=xpix) {x=7*x/8;}
ymaximum = memlimit/x;
y= 1<< ((unsigned int) (floor(log(ymaximum-1)/log(2))+1));
if (7*y/8ypix) {y=5*y/8;}
else
if (3*y/4>=ypix) {y=3*y/4;}
else
if (7*y/8>=ypix) {y=7*y/8;}
xmaximum = memlimit/y;
x= 1<< ((unsigned int) (floor(log(xmaximum-1)/log(2))+1));
if (7*x/8 xsize ? 0 : (xpix-xsize)/2);
//y0cut=(y > ysize ? 0 : (ypix-ysize)/2);
}
else
{
xsize=(foundx < xpix ? foundx : xpix );
ysize=(foundy < ypix ? foundy : ypix );
x=foundx;
y=foundy;
// cut the picture from centre
//x0cut=(x > xsize ? 0 : (xpix-xsize)/8);
//y0cut=(y > ysize ? 0 : (ypix-ysize)/8);
}
}
fprintf(stderr,"zeropadding: x=%d\ty=%d\n",x,y);
#ifndef NDEBUG
fprintf(stderr,"memory: %.2f + %.2f = %.2fMB\n",
xpix*ypix/1024.0/1024.0,
y*x2*sizeof(double)/1024.0/1024.0,
(xpix*ypix+y*x2*sizeof(double))/1024.0/1024.0);
#endif
}
}
win_init(xsize,ysize);
ang=getangle( pix, xpix, ypix,
xsize, ysize,
x0cut, y0cut,
x, y, 1);
win_destroy();
fprintf(stderr,"angle=%f\n", -ang*180/PI);
fprintf(stderr,"pixels dx=%.2f, dy=%.2f\n", sin(ang)*ysize, sin(ang)*xsize);
if ((abs(ang*180/PI)<0.04)&&(abs(sin(ang)*ysize)<2)&&(abs(sin(ang)*xsize)<2))
{
fprintf(stderr,"angle very small, setting angle=0\n");
fprintf(stdout,"0");
}
else
{
fprintf(stdout,"%f", -ang*180/PI);
}
return 0;
}
///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
inline double XY(double *pict, int X, int Y,int y,int x)
{
/*
zero zero
| |
00000 00000000 - zero
00000 11110000
00000 11110000
00000 => 11110000
11111
11111
11111
11111
*/
if (y>0)
if (x>=0)
return(pict[y*X+x]);
else
// the lower part starts from zeroth column too!!!
//(Y-1)== max row index
// the row at the bottom is the first row for x<0
//return(pict[((Y-1)-(y-1))*X-x]); // equivalent is:
return(pict[(Y-y)*X-x]); // y>=1
else
{
//zeroth row is symmetrical
if (y==0) {return(pict[y*X+abs(x)]);}
if (x>=0)
return(pict[-y*X+x]);
else
return(pict[(Y+y)*X-x]);
}
}
static double angle(double *pict, int X, int Y, int Xpix, int Ypix,
double astart, double aend, int anum)
{
int r,X2=(X/2+1);
double d,alpha,x,y,fx,fy,dx1,dx2,dy1,dy2, alpha_max=0;
double s,smax=0;
double w[2][2];
double YY=(Y > X ? Y/2 : X/2) ; //set number of points to bigger axis
double step=(aend-astart)/(anum-1);
#ifndef NDEBUG
fprintf(stderr,"%f\t%f\nstep =%f\n",astart*180/PI,aend*180/PI,step*180/PI);
#endif
alpha=astart;
for(r=0;r X) {sina*=(double) X/(double) Y;}
else {cosa*=(double) Y/(double) X;}
s=0;
for(d=YY;d>3;d-=1) //start from smaller numbers
{
x=d*sina;
y=d*cosa;
fx=floor(x);
fy=floor(y);
dx1=(fx-x)*Xpix/X;
dy1=(fy-y)*Ypix/Y;
dx2=(dx1+1)*Xpix/X;
dy2=(dy1+1)*Ypix/Y;
/*
dx1=(fx-x);
dy1=(fy-y);
dx2=(dx1+1);
dy2=(dy1+1);
*/
w[0][0]=asinc(dx1,X)*asinc(dy1,Y);
w[1][0]=asinc(dx2,X)*asinc(dy1,Y);
w[0][1]=asinc(dx1,X)*asinc(dy2,Y);
w[1][1]=asinc(dx2,X)*asinc(dy2,Y);
/*w[][]
00 10
*
01 11
*/
/*
w[0][0]=1;
w[1][0]=1;
w[0][1]=1;
w[1][1]=1;
*/
s+=(XY(pict,X2,Y,(int) fy, (int) fx)*w[0][0]+
XY(pict,X2,Y,(int) fy+1, (int) fx)*w[0][1]+
XY(pict,X2,Y,(int) fy, (int) fx+1)*w[1][0]+
XY(pict,X2,Y,(int) fy+1, (int) fx+1)*w[1][1]);
// s+=stmp*stmp; //compute power
//s+=XY(pict,X,Y,(int) floor(y+0.5), (int) floor(x+0.5));
}
if (s>smax) { smax=s; alpha_max=alpha;}
#ifdef DEBUG
fprintf(stderr,"%f\t%f\n",alpha/PI*180,s);
#endif
}
#ifdef DEBUG
fprintf(stderr,"%f\t%f\n",alpha_max/PI*180,smax);
#endif
return(alpha_max);
}
inline double asinc(double d, int N)
{
//because N is relatively big we use sinc() instead of asinc()
if (d!=0) {
//cosinus window
//superposition of asinc computed for raised cosinus window!!!
if (d==1) {
return((double) fabs((0.5*PI+
sin(d*PI)/d+
0.5*sin((d+1.0)*PI)/(d+1))/PI)
);
}
if (d==-1) {
return((double) fabs((0.5*sin((d-1.0)*PI)/(d-1)+
sin(d*PI)/d+
0.5*PI)/PI)
);
}
return((double) fabs((0.5*sin((d-1.0)*PI)/(d-1)+
sin(d*PI)/d+
0.5*sin((d+1.0)*PI)/(d+1))/PI)
);
//return sin(d*PI)/(N*sin(d*PI/N)); // real asinc
//return sin(d*PI)/(d*PI); // sinc as an
// approximaton of asinc
}
else return (1);
}
double *xwin,*ywin;
static void win_init(unsigned int Nx, unsigned int Ny)
{
unsigned int r;
double PNx=2*PI/Nx;
double PNy=2*PI/Ny;
double Nx2=Nx/2.0;
double Ny2=Ny/2.0;
double N=(Nypict_max) {pict_max=pict[ry*x2+rx];}
if (pict[ry*x2+rx]0) {z= h;} else {z=0;}
// if (r<20) fprintf(stderr,"%d\n", z);
fwrite(&z,1,1,f);
}
fclose(f);
#ifdef DEBUG
fprintf(stderr,"done\n");
#endif
}
else
{fprintf(stderr,"couldn't open file ./absfft.pnm\n");}
}
#endif
{
int ixmax=imax%(x2/2)
, iymax=imax/(x2/2),
nr=10;
double astart,aend,d;
if (!iymax)
{
nr=31;
ang=0;
astart=-PI/4;
aend=PI/4;
}
else
{
// \pm 2*pixel
if (iymax>y/2)
{ // lower half
iymax=(y-1)-iymax+1;
ang=-atan((float)(ixmax)/(float)iymax*y/x);
astart=-atan((float)(ixmax+1+2)/(float)iymax);
aend=-atan((float)(ixmax-1)/(float)iymax);
}
else
{ //upper half
ang=atan((float)(ixmax)/(float)iymax*y/x);
astart=atan((float)(ixmax-2)/(float)iymax);
aend=atan((float)(ixmax+2)/(float)iymax);
}
}
#ifdef DEBUG
fprintf(stderr,"max %d\t%d\t%f\n", iymax, ixmax,ang*180/PI);
#endif
do
{
ang=angle(pict,x,y,xsize,ysize,astart,aend,nr);
d=(aend-astart)/5;
astart=ang-d*0.66;
aend=ang+d*0.66;
nr=10;
} while(d>0.01*PI/180);
// ang=atan(tan(ang)*((double)y)/((double)x));
if (ang>PI/4) ang-=PI/2;
if (ang<-PI/4) ang+=PI/2;
}
free(pict);
return(ang);
}