彩色图象是多光谱图象的一种特殊情况,对应于人类视觉的三基色即红、绿、蓝三个波段,是对人眼的光谱量化性质的近似。三基色模型是建立图象成像、显示、打印等设备的基础,具有十分重要的作用。在第二章中,已经就彩色做了基本介绍,例如各种彩色空间的表示、定量关系,彩色立方体等直观的图示,参见:(第二章 图象获取、显示、表示与处理)图象显示部分,这里不再重复。
有关彩色的资料参见:efg's Color Page: A. Color Information B. Color Quantization C. Palettes 。
彩色空间有多种不同的表示形式,常用的有红、绿、蓝(RGB)三基色模型,色调、饱和度、亮度(HSI)模型。参见:(第二章 图象获取、显示、表示与处理)图象显示部分。
美国国家电视系统委员会(NTSC)定义了用光亮度和色度传送信号的格式YIQ,其中Y代表亮度信息,I、Q为色度值[1];欧洲定义了相交替格式(Phase Alternating Line (PAL)),使用YUV格式。YUV格式与YIQ格式类似,差别仅在于空间上多一个330的旋转。与YUV彩色空间具有数字等价性的YCbCr彩色空间是以演播室质量标准为目标的CCIR601编码方案中采用的彩色表示模型,Cb与U分量对应,而Cr与V分量对应。在该编码方案中,亮度信号与色度信号Cb、Cr的采样比率为4:2:2,这是因为人眼对色度信号的变化没有对亮度信号的变化来得敏感。
青、品红、黄(CMY)彩色模型[3]是彩色图象印刷行业使用的彩色空间,在彩色立方体中它们是红、绿、蓝的补色,称为减色基,而红、绿、蓝称为加色基。在CMY模型中,颜色是从白光中减去一定成分得到的。CMY坐标可以从RGB模型中得到:
由于在印刷时CMY模型不可能产生真正的黑色,因此在印刷业中实际上使用的是CMYK彩色模型,K为第四种颜色,表示黑色:从CMY 到CMYK的转换:
此外,还有国际照明委员会(CIE)提出的基于色度学的彩色模型[4]。
本节的公式一般是指在归一化坐标下的转换公式,即亮度数值最大为1。
3.1 RGB与HSI
对于彩色图象分割而言,有时需要将RGB变换为HSI坐标,以便反映人类观察彩色的方式,转换公式如下:
RGB与HSI的转换关系:图形显示程序efg's HSV Lab Report
反过来,由HSI变换为RGB坐标的公式如下:
3.2 RGB与YIQ
3.3 RGB与YUV(YCbCr)
特别地对于RGB(8位表示的RGB)与YCbCr(256级别)之间的转换如下:
YCbCr (256 级别) 可以从8位 RGB 直接计算:
Y = 0.299 R + 0.587 G + 0.114 B
Cb = - 0.1687 R - 0.3313 G + 0.5 B + 128
Cr = 0.5 R - 0.4187 G - 0.0813 B + 128
YCbCr to RGB Conversion
反过来,RGB 也可以直接从YCbCr (256级别) 计算:
R = Y + 1.402 (Cr-128)
G = Y - 0.34414 (Cb-128) - 0.71414 (Cr-128)
B = Y + 1.772 (Cb-128)
彩色图象的分割除了可以采用灰度图象的分割方法外,还可以利用色彩的聚类性质。通常物体的色度和饱和度由构成物体的材料所具有的光线吸收和反射特性决定,而亮度明显地受光照和视角的影响。因此根据色度和饱和度分割图象比较可靠。在实际问题中,需要对一类物体的颜色进行统计训练,考察其颜色的直方图,根据颜色直方图选择合适的范围可以比较快速地实现图象分割。
图象的几何变换
这一章我们将介绍图象的几何变换,包括图象的平移、旋转、镜象变换、转置、放缩等。如果你熟悉矩阵运算,你将发现,实现这些变换是非常容易的。
2.1 平移
平移(translation)变换大概是几何变换中最简单的一种了。
如图2.1所示,初始坐标为(x0,y0)的点经过平移(tx,ty)(以向右,向下为正方向)后,坐标变为(x1,y1)。这两点之间的关系是x1=x0+tx ,y1=y0+ty。
图2.1 平移的示意图
以矩阵的形式表示为
(2.1)
我们更关心的是它的逆变换:
(2.2)
这是因为:我们想知道的是平移后的图象中每个象素的颜色。例如我们想知道,新图中左上角点的RGB值是多少?很显然,该点是原图的某点经过平移后得到的,这两点的颜色肯定是一样的,所以只要知道了原图那点的RGB值即可。那么到底新图中的左上角点对应原图中的哪一点呢?将左上角点的坐标(0,0)入公式(2.2),得到x0=-tx ,y0=-ty;所以新图中的(0,0)点的颜色和原图中(-tx , -ty)的一样。
这样就存在一个问题:如果新图中有一点(x1,y1),按照公式(2.2)得到的(x0,y0)不在原图中该怎么办?通常的做法是,把该点的RGB值统一设成(0,0,0)或者(255,255,255)。
另一个问题是:平移后的图象是否要放大?一种做法是不放大,移出的部分被截断。例如,图2.2为原图,图2.3为移动后的图。这种处理,文件大小不会改变。
图2.2 移动前的图
图2.3 移动后的图
还有一种做法是:将图象放大,使得能够显示下所有部分,如图2.4所示。
图2.4 移动后图象被放大
这种处理,文件大小要改变。设原图的宽和高分别是w1,h1则新图的宽和高变为w1+|tx|和h1+|ty|,加绝对值符号是因为tx, ty有可能为负(即向左,向上移动)。
下面的函数Translation采用的是第一种做法,即移出的部分被截断。在给出源代码之前,先说明一个问题。
如果你用过Photoshop,Corel PhotoPaint等图象处理软件,可能听说过“灰度图”(grayscale)这个词。灰度图是指只含亮度信息,不含色彩信息的图象,就象我们平时看到的黑白照片:亮度由暗到明,变化是连续的。因此,要表示灰度图,就需要把亮度值进行量化。通常划分成0到255共256个级别,其中0最暗(全黑),255最亮(全白)。.bmp格式的文件中,并没有灰度图这个概念,但是,我们可以很容易在.bmp文件中表示灰度图。方法是用256色的调色板,只不过这个调色板有点特殊,每一项的RGB值都是相同的。也就是说RGB值从(0,0,0),(1,1,1)一直到(255,255,255)。(0,0,0)是全黑色,(255,255,255)是全白色,中间的是灰色。这样,灰度图就可以用256色图来表示了。为什么会这样呢?难道是一种巧合?其实并不是。
在表示颜色的方法中,除了RGB外,还有一种叫YUV的表示方法,应用也很多。电视信号中用的就是一种类似于YUV的颜色表示方法。
在这种表示方法中,Y分量的物理含义就是亮度,U和V分量代表了色差信号(你不必了解什么是色差,只要知道有这么一个概念就可以了)。使用这种表示方法有很多好处,最主要的有两点:
(1) 因为Y代表了亮度,所以Y分量包含了灰度图的所有信息,只用Y分量就能完全能够表示出一幅灰度图来。当同时考虑U,V分量时,就能够表示出彩色信息来。这样,用同一种表示方法可以很方便的在灰度和彩色图之间切换,而RGB表示方法就做不到这一点了。
(2) 人眼对于亮度信号非常敏感,而对色差信号的敏感程度相对较弱。也就是说,图象的主要信息包含在Y分量中。这就提示我们:如果在对YUV信号进行量化时,可以“偏心”一点,让Y的量化级别多一些(谁让它重要呢?)而让UV的量化级别少一些,就可以实现图象信息的压缩。这一点将在第9章介绍图象压缩时仔细研究,这里就不深入讨论了。而RGB的表示方法就做不到这一点,因为RGB三个分量同等重要,缺了谁也不行。YUV和RGB之间有着如下的对应关系
(2.3)
(2.4)
当RGB三个分量的大小一样时,假设都是a,代入公式(2.3),得到Y=a,U=0,V=0 。你现在该明白我前面所说不是巧合的原因了吧。
使用灰度图有一个好处,那就是方便。首先RGB的值都一样;其次,图象数据即调色板索引值,也就是实际的RGB值,也就是亮度值;另外,因为是256色调色板,所以图象数据中一个字节代表一个象素,很整齐。如果是2色图或16色图,还要拼凑字节,很麻烦。如果是彩色的256色图,由于图象处理后有可能会产生不属于这256种颜色的新颜色,就更麻烦了;这一点,今后你就会有深刻体会的。所以,做图象处理时,一般采用灰度图。为了将重点放在算法本身上,今后给出的程序如不做特殊说明,都是针对256级灰度图的。其它颜色的情况,你可以自己想一想,把算法补全。
如果想得到一幅灰度图,可以使用Sea或者PhotoShop等软件提供的颜色转换功能将彩色图转换成灰度图。
好了,言归正传,下面给出Translation的源代码。算法的思想是先将所有区域填成白色,然后找平移后显示区域的左上角点(x0,y0) 和右下角点(x1,y1) ,分几种情况进行处理。
先看x方向(width指图象的宽度)
(1) tx≤-width:很显然,图象完全移出了屏幕,不用做任何处理;
(2) -width<tx≤0:如图2.5所示。容易看出,图象区域的x范围从0到width-|tx|,对应原图的范围从|tx|到width;
图2.5 tx≤0,ty≤0的情况
(3) 0< tx <width:如图2.6所示。容易看出,图象区域的x范围从tx 到width,对应原图的范围从0到width - tx ;
图2.6 0< tx<width,0<ty<height的情况
(4) tx ≥width:很显然,图象完全移出了屏幕,不用做任何处理。
y方向是对应的(height表示图象的高度):
(1) ty≤-height,图象完全移出了屏幕,不用做任何处理;
(2) -height<ty≤0,图象区域的y范围从0到height-|ty|,对应原图的范围从|ty|到height;
(3) 0<ty<height ,图象区域的y范围从ty到height,对应原图的范围从0到height-ty;
(4) ty≥height,图象完全移出了屏幕,不用做任何处理。
这种做法利用了位图存储的连续性,即同一行的象素在内存中是相邻的。利用memcpy函数,从(x0,y0)点开始,一次可以拷贝一整行(宽度为x1-x0),然后将内存指针移到(x0,y0+1)处,拷贝下一行。这样拷贝(y1-y0)行就完成了全部操作,避免了一个一个象素的计算,提高了效率。Translation的源代码如下:
int xOffset=0,yOffset=0;
BOOL Translation(HWND hWnd)
{
DLGPROC dlgInputBox = NULL;
DWORD OffBits,BufSize;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
int SrcX0,SrcY0,SrcX1,SrcY1;
int DstX0,DstY0,DstX1,DstY1;
int RectWidth,RectHeight;
BOOL xVisible,yVisible;
HDC hDc;
HFILE hf;
int i;
//出现对话框,输入x偏移量xOffset,和y偏移量yOffset
dlgInputBox = (DLGPROC) MakeProcInstance ( (FARPROC)InputBox,ghInst );
DialogBox (ghInst, "INPUTBOX", hWnd, dlgInputBox);
FreeProcInstance ( (FARPROC) dlgInputBox );
//OffBits为BITMAPINFOHEADER结构长度加调色板的大小
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
BufSize=OffBits+bi.biHeight*LineBytes;//要开的缓冲区的大小
//为新产生的位图分配缓冲区内存
if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|
MB_ICONEXCLAMATION);
return FALSE; //失败,返回
}
//lpImgData为指向原来位图数据的指针
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);
//lpTempImgData为指向新产生位图数据的指针
lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
lpPtr=(char *)lpImgData;
lpTempPtr=(char *)lpTempImgData;
//将新的缓冲区中的每个字节都填成255,这样以后未处理的象素就是白色
memset(lpTempPtr,(BYTE)255,BufSize);
//两幅图之间的头信息,包括调色板都是相同的,所以直接拷贝头和调色板
memcpy(lpTempPtr,lpPtr,OffBits);
//xVisible为FALSE时,表示x方向已经移出了可显示的范围
xVisible=TRUE;
if( xOffset<= -bi.biWidth )
xVisible=FALSE;
else if( xOffset<=0){
DstX0=0; //表示移动后,有图区域的左上角点的x坐标
DstX1=bi.biWidth+xOffset; //表示移动后,有图区域的右下角点的x坐标
}
else if ( xOffset<bi.biWidth){
DstX0=xOffset;
DstX1=bi.biWidth;
}
else
xVisible=FALSE;
SrcX0=DstX0-xOffset; //对应DstX0在原图中的x坐标
SrcX1=DstX1-xOffset; //对应DstX1在原图中的x坐标
RectWidth=DstX1-DstX0; //有图区域的宽度
//yVisible为FALSE时,表示y方向已经移出了可显示的范围
yVisible=TRUE;
if( yOffset<= -bi.biHeight )
yVisible=FALSE;
else if( yOffset<=0){
DstY0=0; //表示移动后,有图区域的左上角点的y坐标
DstY1=bi.biHeight+yOffset; //表示移动后,有图区域的右下角点的y坐标
}
else if ( yOffset<bi.biHeight){
DstY0=yOffset;
DstY1=bi.biHeight;
}
else
yVisible=FALSE;
SrcY0=DstY0-yOffset; //对应DstY0在原图中的y坐标
SrcY1=DstY1-yOffset; //对应DstY1在原图中的y坐标
RectHeight=DstY1-DstY0; //有图区域的高度
if( xVisible && yVisible){ //x,y方向都没有完全移出可显示的范围
for(i=0;i<RectHeight;i++){ //拷贝每一行
//lpPtr指向要拷贝的那一行的最左边的象素对应在原图中的位
//置。特别要注意的是,由于.bmp是上下颠倒的,偏移是
//(BufSize-LineBytes-(i+SrcY0)*LineBytes)+SrcX0,而不是
//(i+SrcY0)*LineBytes)+SrcX0,你试着举个例子就明白了。
lpPtr=(char*)lpImgData+(BufSize-LineBytes-
(i+SrcY0)*LineBytes)+SrcX0;
//lpTempPtr指向要拷贝的那一行的最左边的象素对应在新图中//的位置。同样要注意上面//的问题。
lpTempPtr=(char*)lpTempImgData+
(BufSize-LineBytes-(i+DstY0)*LineBytes)+DstX0;
//拷贝一行(宽度为RectWidth)
memcpy(lpTempPtr,lpPtr,RectWidth);
}
}
hDc=GetDC(hWnd);
if(hBitmap!=NULL)
DeleteObject(hBitmap); //释放原来的位图句柄
//产生新的位图
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER) +
NumColors*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData,
DIB_RGB_COLORS);
//将平移后的图象存成文件
hf=_lcreat("c:\\translation.bmp",0);
_lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,BufSize);
_lclose(hf);
//释放资源和内存
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}
2.2 旋转
旋转(rotation)有一个绕着什么转的问题,通常的做法是以图象的中心为圆心旋转,举个例子,图2.7旋转30度(顺时针方向)后如图2.8所示:
图2.7 旋转前的图 |
图2.8 旋转后的图 |
可以看出,旋转后图象变大了。另一种做法是不让图象变大,转出的部分被裁剪掉。如图2.9所示。
我们采用第一种做法,首先给出变换矩阵。在我们熟悉的坐标系中,将一个点顺时针旋转a角后的坐标变换公式,如图2.10所示,r为该点到原点的距离,在旋转过程中,r保持不变;b为r与x轴之间的夹角。
图2.9 旋转后保持原图大小, 转出的部分被裁掉 |
图2.10 旋转示意图 |
旋转前:x0=rcosb;y0=rsinb
旋转a角度后:
x1=rcos(b-a)=rcosbcosa+rsinbsina=x0cosa+y0sina;
y1=rsin(b-a)=rsinbcosa-rcosbsina=-x0sina+y0cosa;
以矩阵的形式表示:
(2.5)
上面的公式中,坐标系xoy是以图象的中心为原点,向右为x轴正方向,向上为y轴正方向。它和以图象左上角点为原点o’,向右为x’轴正方向,向下为y’轴正方向的坐标系x’o’y’之间的转换关系如何呢?如图2.11所示。
图2.11 两种坐标系间的转换关系
设图象的宽为w,高为h,容易得到:
(2.6)
逆变换为:
(2.7)
有了上面的公式,我们可以把变换分成三步:
1.将坐标系o’变成o;
2.将该点顺时针旋转a角;
3.将坐标系o变回o’,这样,我们就得到了变换矩阵,是上面三个矩阵的级联。
(2.8)
要注意的是,因为新图变大,所以上面公式中出现了wold,hold,wnew,hnew,它们分别表示原图(old)和新图(new)的宽、高。我们从图2.8中容易看出:wnew=max(|x4-x1|,|x3-x2|) ;hnew=max(|y4-y1|,|y3-y2|)。
(2.8)的逆变换为
(2.9)
这样,对于新图中的每一点,我们就可以根据公式(2.9)求出对应原图中的点,得到它的灰度。如果超出原图范围,则填成白色。要注意的是,由于有浮点运算,计算出来点的坐标可能不是整数,采用取整处理,即找最接近的点,这样会带来一些误差(图象可能会出现锯齿)。更精确的方法是采用插值,将在图象缩放时介绍。
源程序如下:
#define PI 3.1415926535
#define RADIAN(angle) ((angle)*PI/180.0) //角度到弧度转化的宏
BOOL Rotation(HWND hWnd)
{
DLGPROC dlgInputBox = NULL;
DWORD OffBits,SrcBufSize,DstBufSize,DstLineBytes;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
float SrcX1,SrcY1,SrcX2,SrcY2;
float SrcX3,SrcY3,SrcX4,SrcY4;
float DstX1,DstY1,DstX2,DstY2;
float DstX3,DstY3,DstX4,DstY4;
DWORD Wold,Hold,Wnew,Hnew;
HDC hDc;
HFILE hf;
DWORD x0,y0,x1,y1;
float cosa,sina; //cos(a),sin(a);
float num1,num2;
BITMAPFILEHEADER DstBf;
BITMAPINFOHEADER DstBi;
//出现对话框,输入旋转角度(顺时针方向)
dlgInputBox = (DLGPROC) MakeProcInstance ( (FARPROC)InputBox,ghInst );
DialogBox (ghInst, "INPUTBOX", hWnd, dlgInputBox);
FreeProcInstance ( (FARPROC) dlgInputBox );
//角度到弧度的转化
RotateAngle=(float)RADIAN(RotateAngle);
cosa=(float)cos((double)RotateAngle);
sina=(float)sin((double)RotateAngle);
//原图的宽度和高度
Wold=bi.biWidth;
Hold=bi.biHeight;
//原图的四个角的坐标
SrcX1=(float)(-0.5*Wold);
SrcY1=(float)(0.5*Hold);
SrcX2=(float)(0.5*Wold);
SrcY2=(float)(0.5*Hold);
SrcX3=(float)(-0.5*Wold);
SrcY3=(float)(-0.5*Hold);
SrcX4=(float)(0.5*Wold);
SrcY4=(float)(-0.5*Hold);
//新图四个角的坐标
DstX1=cosa*SrcX1+sina*SrcY1;
DstY1=-sina*SrcX1+cosa*SrcY1;
DstX2=cosa*SrcX2+sina*SrcY2;
DstY2=-sina*SrcX2+cosa*SrcY2;
DstX3=cosa*SrcX3+sina*SrcY3;
DstY3=-sina*SrcX3+cosa*SrcY3;
DstX4=cosa*SrcX4+sina*SrcY4;
DstY4=-sina*SrcX4+cosa*SrcY4;
//计算新图的宽度,高度
Wnew = (DWORD)(max(fabs(DstX4-DstX1), fabs(DstX3-DstX2))+0.5);
Hnew = (DWORD)(max(fabs(DstY4-DstY1), fabs(DstY3-DstY2))+0.5);
//计算矩阵(2.9)中的两个常数,这样不用以后每次都计算了
num1=(float)( -0.5*Wnew*cosa-0.5*Hnew*sina+0.5*Wold);
num2=(float)(0.5*Wnew*sina-0.5*Hnew*cosa+0.5*Hold);
//OffBits为BITMAPINFOHEADER结构长度加调色板的大小
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
SrcBufSize=OffBits+bi.biHeight*LineBytes;
//显示时,采用新图的宽度和高度,
ImgWidth=Wnew;
ImgHeight=Hnew;
//新图每行占用的字节
DstLineBytes=(DWORD)WIDTHBYTES(Wnew*bi.biBitCount);
DstBufSize=(DWORD)(sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD)+
(DWORD)DstLineBytes*Hnew); //要开的缓冲区的大小
//为新产生的位图分配缓冲区内存
if((hTempImgData=LocalAlloc(LHND,DstBufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",
MB_OK|MB_ICONEXCLAMATION);
return FALSE; //失败,返回
}
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);
//lpTempImgData为指向新产生位图数据的指针
lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
lpPtr=(char *)lpImgData;
lpTempPtr=(char *)lpTempImgData;
//将新的缓冲区中的每个字节都填成255,这样以后未处理的象素就是白色
memset(lpTempPtr,(BYTE)255,DstBufSize);
//拷贝头和调色板信息
memcpy(lpTempPtr,lpPtr,OffBits);
//得到新的BITMAPFILEDER和BITMAPINFOHERDER
memcpy((char *)&DstBf,(char *)&bf,sizeof(BITMAPFILEHEADER));
memcpy((char *)&DstBi,(char *)&bi,sizeof(BITMAPINFOHEADER));
//做一些必要的改变,这一点特别要注意
DstBf.bfSize=DstBufSize+sizeof(BITMAPFILEHEADER);
DstBi.biWidth=Wnew;
DstBi.biHeight=Hnew;
//用新的BITMAPINFOHERDER覆盖原来的那个
memcpy(lpTempPtr,(char *)&DstBi,sizeof(BITMAPINFOHEADER));
for(y1=0;y1<Hnew;y1++)
for(x1=0;x1<Wnew;x1++){
//x0,y0为对应的原图上的坐标
x0= (DWORD)(x1*cosa+y1*sina+num1);
y0= (DWORD)(-1.0f*x1*sina+y1*cosa+num2);
if( (x0>=0) && (x0<Wold) && (y0>=0) && (y0<Hold))
//在原图范围内
{
lpPtr=(char*)lpImgData+
(SrcBufSize-LineBytes-y0*LineBytes)+x0;
lpTempPtr=(char*)lpTempImgData+
(DstBufSize-DstLineBytes-y1*DstLineBytes)+x1;
*lpTempPtr=*lpPtr; //进行象素的复制
}
}
hDc=GetDC(hWnd);
if(hBitmap!=NULL)
DeleteObject(hBitmap); //释放原来的位图句柄
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER) +
NumColors*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData,
DIB_RGB_COLORS);
//将旋转后的图象存成文件
hf=_lcreat("c:\\rotation.bmp",0);
_lwrite(hf,(LPSTR)&DstBf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,DstBufSize);
_lclose(hf);
//释放资源和内存
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}
程序运行时的画面如图2.12所示
图2.12 旋转
2.3 镜象
镜象(mirror)分水平镜象和垂直镜象两种。图2.2的水平镜象和垂直镜象分别如图2.13和图2.14所示
图2.13 图2.2的水平镜象
图2.14 图2.2的垂直镜象
镜象的变换矩阵很简单。设原图宽为w,高为h,变换后,图的宽和高不变。
水平镜象的变化矩阵为:
(2.10)
垂直镜象的变化矩阵为:
(2.11)
镜象变换的源代码如下,因为和平移的那段程序很类似,程序中的注释就简单一些。
BOOL Mirror(HWND hWnd,BOOL XDirection)
//Xdirection为TRUE时表示水平镜象,为FALSE时表示垂直镜象变换
{
DWORD OffBits,BufSize;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
HDC hDc;
HFILE hf;
LONG x0,y0,x1,y1;
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
BufSize=OffBits+bi.biHeight*LineBytes;
if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|
MB_ICONEXCLAMATION);
return FALSE;
}
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData); lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
lpPtr=(char *)lpImgData;
lpTempPtr=(char *)lpTempImgData;
memset(lpTempPtr,(BYTE)255,BufSize);
memcpy(lpTempPtr,lpPtr,OffBits);
if( XDirection){ //水平镜象
for(y1=0;y1<bi.biHeight;y1++)
for(x1=0;x1<bi.biWidth;x1++){
x0=bi.biWidth-1-x1; //因为x坐标是从0到bi.biWidth-1
y0=y1;
lpPtr=(char *)lpImgData+(BufSize-LineBytes-y0*LineBytes)+x0;
lpTempPtr=(char *)lpTempImgData+
(BufSize-LineBytes-y1*LineBytes)+x1;
*lpTempPtr=*lpPtr;
}
}
else{ //垂直镜象
for(y1=0;y1<bi.biHeight;y1++)
for(x1=0;x1<bi.biWidth;x1++){
x0=x1;
y0=bi.biHeight-1-y1;
lpPtr=(char *)lpImgData+(BufSize-LineBytes-y0*LineBytes)+x0;
lpTempPtr=(char *)lpTempImgData+
(BufSize-LineBytes-y1*LineBytes)+x1;
*lpTempPtr=*lpPtr;
}
}
hDc=GetDC(hWnd);
if(hBitmap!=NULL)
DeleteObject(hBitmap);
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData, DIB_RGB_COLORS);
if( XDirection)
hf=_lcreat("c:\\mirrorx.bmp",0);
else
hf=_lcreat("c:\\mirrory.bmp",0);
_lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,BufSize);
_lclose(hf);
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}
2.4 转置
转置(transpose)是指将x,y坐标对换,图2.2的转置如图2.15所示。
图2.15 图2.2的转置
要注意的是,转置和旋转900是有区别的,不信你可以试试:怎么旋转,图2.2也转不出图2.15来。另外,转置后图的宽高对换了。转置的变换矩阵很简单:
(2.12)
镜象变换的源代码如下,因为和旋转的那段程序很类似,程序中的注释就简单一些:
BOOL Transpose(HWND hWnd)
{
DWORD OffBits,SrcBufSize,DstBufSize,DstLineBytes;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
DWORD Wnew,Hnew;
HDC hDc;
HFILE hf;
DWORD x0,y0,x1,y1;
BITMAPFILEHEADER DstBf;
BITMAPINFOHEADER DstBi;
//新图的宽度和高度
Wnew = (DWORD)bi.biHeight;
Hnew = (DWORD)bi.biWidth;
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
SrcBufSize=OffBits+bi.biHeight*LineBytes;
//显示时,采用新图的宽度和高度,
ImgWidth=Wnew;
ImgHeight=Hnew;
DstLineBytes=(DWORD)WIDTHBYTES(Wnew*bi.biBitCount);
DstBufSize=(DWORD)(sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD)+
(DWORD)DstLineBytes*Hnew);
if((hTempImgData=LocalAlloc(LHND,DstBufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|
MB_ICONEXCLAMATION);
return FALSE;
}
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);
lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
lpPtr=(char *)lpImgData;
lpTempPtr=(char *)lpTempImgData;
memset(lpTempPtr,(BYTE)255,DstBufSize);
memcpy(lpTempPtr,lpPtr,OffBits);
//头信息中做一些必要的改变,这一点非常重要
memcpy((char *)&DstBf,(char *)&bf,sizeof(BITMAPFILEHEADER));
memcpy((char *)&DstBi,(char *)&bi,sizeof(BITMAPINFOHEADER));
DstBf.bfSize=DstBufSize+sizeof(BITMAPFILEHEADER);
DstBi.biWidth=Wnew;
DstBi.biHeight=Hnew;
memcpy(lpTempPtr,(char *)&DstBi,sizeof(BITMAPINFOHEADER));
for(y1=0;y1<Hnew;y1++)
for(x1=0;x1<Wnew;x1++){
x0= y1;
y0= x1;
lpPtr=(char *)lpImgData+(SrcBufSize-LineBytes-y0*LineBytes)+x0;
lpTempPtr=(char *)lpTempImgData+
(DstBufSize-DstLineBytes-y1*DstLineBytes)+x1;
*lpTempPtr=*lpPtr;
}
hDc=GetDC(hWnd);
if(hBitmap!=NULL)
DeleteObject(hBitmap);
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData, DIB_RGB_COLORS);
hf=_lcreat("c:\\transpose.bmp",0);
_lwrite(hf,(LPSTR)&DstBf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,DstBufSize);
_lclose(hf);
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}
2.5 缩放
假设放大因子为ratio,(为了避免新图过大或过小,我们在程序中限制0.25≤ratio≤4),缩放(zoom)的变换矩阵很简单:
(2.13)
缩放变换的源代码如下,因为和转置的那段程序很类似,程序中的注释就简单一些。
float ZoomRatio=0.25f; //缩放比例,初始化为0.25
BOOL Zoom(HWND hWnd)
{
DLGPROC dlgInputBox = NULL;
DWORD OffBits,SrcBufSize,DstBufSize,DstLineBytes;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
DWORD Wold,Hold,Wnew,Hnew;
HDC hDc;
HFILE hf;
DWORD x0,y0,x1,y1;
float num1;
BITMAPFILEHEADER DstBf;
BITMAPINFOHEADER DstBi;
//出现对话框,输入缩放比例
dlgInputBox = (DLGPROC) MakeProcInstance ( (FARPROC)InputBox, ghInst );
DialogBox (ghInst, "INPUTBOX", hWnd, dlgInputBox);
FreeProcInstance ( (FARPROC) dlgInputBox );
num1=(float)(1.0/ZoomRatio);
//原图宽度和高度
Wold=bi.biWidth;
Hold=bi.biHeight;
//新图宽度和高度
Wnew = (DWORD)(Wold*ZoomRatio+0.5);
Hnew = (DWORD)(Hold*ZoomRatio+0.5);
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
SrcBufSize=OffBits+bi.biHeight*LineBytes;
ImgWidth=Wnew;
ImgHeight=Hnew;
DstLineBytes=(DWORD)WIDTHBYTES(Wnew*bi.biBitCount);
DstBufSize=(DWORD)(sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD)+
(DWORD)DstLineBytes*Hnew);
if((hTempImgData=LocalAlloc(LHND,DstBufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|
MB_ICONEXCLAMATION);
return FALSE;
}
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);
lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
lpPtr=(char *)lpImgData;
lpTempPtr=(char *)lpTempImgData;
memset(lpTempPtr,(BYTE)255,DstBufSize);
memcpy(lpTempPtr,lpPtr,OffBits);
//头信息中做一些必要的改变,这一点非常重要
memcpy((char *)&DstBf,(char *)&bf,sizeof(BITMAPFILEHEADER));
memcpy((char *)&DstBi,(char *)&bi,sizeof(BITMAPINFOHEADER));
DstBf.bfSize=DstBufSize+sizeof(BITMAPFILEHEADER);
DstBi.biWidth=Wnew;
DstBi.biHeight=Hnew;
memcpy(lpTempPtr,(char *)&DstBi,sizeof(BITMAPINFOHEADER));
for(y1=0;y1<Hnew;y1++)
for(x1=0;x1<Wnew;x1++){
x0= (DWORD)(x1*num1);
y0= (DWORD)(y1*num1);
if( (x0>=0) && (x0<Wold) && (y0>=0) && (y0<Hold))
{
lpPtr=(char*)lpImgData+
(SrcBufSize-LineBytes-y0*LineBytes)+x0;
lpTempPtr=(char *)lpTempImgData+
(DstBufSize-DstLineBytes-y1*DstLineBytes)+x1;
*lpTempPtr=*lpPtr;
}
}
hDc=GetDC(hWnd);
if(hBitmap!=NULL)
DeleteObject(hBitmap);
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER)+
NumColors*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData, DIB_RGB_COLORS);
hf=_lcreat("c:\\zoom.bmp",0);
_lwrite(hf,(LPSTR)&DstBf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,DstBufSize);
_lclose(hf);
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}
由于放大图象时产生了新的象素,以及浮点数的操作,得到的坐标可能并不是整数,这一点我们在介绍旋转时就提到了。我们采用的做法是找与之最临近的点。实际上,更精确的做法是采用插值(interpolation),即利用邻域的象素来估计新的象素值。其实我们前面的做法也是一种插值,称为最邻近插值(Nearest Neighbour Interpolation)。下面先介绍线形插值(Linear Interpolation)。
线形插值使用原图中两个值来构造所求坐标处的值。举一个一维的例子。如图2.16所示,如果已经知道了两点x0,x2处的函数值f(x0),f(x2),现在要求x1处的函数值f(x1)。我们假设函数是线形的,利用几何知识可以知道
f(x1)=(f(x2)-f(x0))(x1-x0)/(x2-x0)+f(x0)
(2.13)
在图象处理中需要将线形插值扩展到二维的情况,即采用双线形插值(Bilinear Intrepolation),图2.17为双线形插值的示意图。
图2.16 线形插值的示意图 |
图.217 双线形插值的示意图 |
已知a、b、c、d四点的灰度,要求e点的灰度,可以先在水平方向上由a,b线形插值求出g、c、d线形插值求出f,然后在垂直方向上由g,f线形插值求出e。
线形插值基于这样的假设:原图的灰度在两个象素之间是线形变化的。一般情况下,这种插值的效果还不错。更精确的方法是采用曲线插值(Curvilinear Interpolation),即认为象素之间的灰度变化规律符合某种曲线,但这种处理的计算量是很大的。
关于插值,我们就介绍到这里,有兴趣的读者可以参考“数值分析”方面的书籍。
RGB to YCbCr
void rgb2ycbcr(int r, int g, int b, int[] ycbcr) {
int y = (int)( 0.299 * r + 0.587 * g + 0.114 * b);
int cb = (int)(-0.16874 * r - 0.33126 * g + 0.50000 * b);
int cr = (int)( 0.50000 * r - 0.41869 * g - 0.08131 * b);
ycbcr[0] = y;
ycbcr[1] = cb;
ycbcr[2] = cr;
}
RGB to YUV
void rgb2yuv(int r,int g, int b, int[] yuv,) {
int y = (int)(0.299 * r + 0.587 * g + 0.114 * b);
int u = (int)((b - y) * 0.492f);
int v = (int)((r - y) * 0.877f);
yuv[0]= y;
yuv[1]= u;
yuv[2]= v;
}RGB to HSB
void rgb2hsb(int r, int g, int b, int[] hsb) {
float [] hsbvals = new float[3];
Color.RGBtoHSB(r, g, b, hsbvals);
}
RGB to HMMD
void rgb2hmmd(int r, int g, int b, int[] hmmd) {
float max = (int)Math.max(Math.max(r,g), Math.max(g,b));
float min = (int)Math.min(Math.min(r,g), Math.min(g,b));
float diff = (max - min);
float sum = (float) ((max + min)/2.);
float hue = 0;
if (diff == 0)
hue = 0;
else if (r == max && (g - b) > 0)
hue = 60*(g-b)/(max-min);
else if (r == max && (g - b) <= 0)
hue = 60*(g-b)/(max-min) + 360;
else if (g == max)
hue = (float) (60*(2.+(b-r)/(max-min)));
else if (b == max)
hue = (float) (60*(4.+(r-g)/(max-min)));
hmmd[0] = (int)(hue); hmmd[1] = (int)(max); hmmd[2] = (int)(min); hmmd[3] = (int)(diff);
}
RGB to HSL
private void rgb2hsl(int r, int g, int b, int hsl[]) {
float var_R = ( r / 255f );
float var_G = ( g / 255f );
float var_B = ( b / 255f );
float var_Min; //Min. value of RGB
float var_Max; //Max. value of RGB
float del_Max; //Delta RGB value
if (var_R > var_G)
{ var_Min = var_G; var_Max = var_R; }
else
{ var_Min = var_R; var_Max = var_G; }
if (var_B > var_Max) var_Max = var_B;
if (var_B < var_Min) var_Min = var_B;
del_Max = var_Max - var_Min;
float H = 0, S, L;
L = ( var_Max + var_Min ) / 2f;
if ( del_Max == 0 ) { H = 0; S = 0; } // gray
else { //Chroma
if ( L < 0.5 )
S = del_Max / ( var_Max + var_Min );
else
S = del_Max / ( 2 - var_Max - var_Min );
float del_R = ( ( ( var_Max - var_R ) / 6f ) + ( del_Max / 2f ) ) / del_Max;
float del_G = ( ( ( var_Max - var_G ) / 6f ) + ( del_Max / 2f ) ) / del_Max;
float del_B = ( ( ( var_Max - var_B ) / 6f ) + ( del_Max / 2f ) ) / del_Max;
if ( var_R == var_Max )
H = del_B - del_G;
else if ( var_G == var_Max )
H = ( 1 / 3f ) + del_R - del_B;
else if ( var_B == var_Max )
H = ( 2 / 3f ) + del_G - del_R;
if ( H < 0 ) H += 1;
if ( H > 1 ) H -= 1;
}
hsl[0] = (int)(360*H);
hsl[1] = (int)(S*100);
hsl[2] = (int)(L*100);
}
RGB to HSV
private void rgb2hsv(int r, int g, int b, int hsv[]) {
int min; //Min. value of RGB
int max; //Max. value of RGB
int delMax; //Delta RGB value
if (r > g) { min = g; max = r; }
else { min = r; max = g; }
if (b > max) max = b;
if (b < min) min = b;
delMax = max - min;
float H = 0, S;
float V = max;
if ( delMax == 0 ) { H = 0; S = 0; }
else {
S = delMax/255f;
if ( r == max )
H = ( (g - b)/(float)delMax)*60;
else if ( g == max )
H = ( 2 + (b - r)/(float)delMax)*60;
else if ( b == max )
H = ( 4 + (r - g)/(float)delMax)*60;
}
hsv[0] = (int)(H);
hsv[1] = (int)(S*100);
hsv[2] = (int)(V*100);
}
RGB to xyY
public void rgb2xyY(int R, int G, int B, int []xyy) {
//http://www.brucelindbloom.com
float rf, gf, bf;
float r, g, b, X, Y, Z;
// RGB to XYZ
r = R/255.f; //R 0..1
g = G/255.f; //G 0..1
b = B/255.f; //B 0..1
if (r <= 0.04045)
r = r/12;
else
r = (float) Math.pow((r+0.055)/1.055,2.4);
if (g <= 0.04045)
g = g/12;
else
g = (float) Math.pow((g+0.055)/1.055,2.4);
if (b <= 0.04045)
b = b/12;
else
b = (float) Math.pow((b+0.055)/1.055,2.4);
X = 0.436052025f*r + 0.385081593f*g + 0.143087414f *b;
Y = 0.222491598f*r + 0.71688606f *g + 0.060621486f *b;
Z = 0.013929122f*r + 0.097097002f*g + 0.71418547f *b;
float x;
float y;
float sum = X + Y + Z;
if (sum != 0) {
x = X / sum;
y = Y / sum;
}
else {
float Xr = 0.964221f; // reference white
float Yr = 1.0f;
float Zr = 0.825211f;
x = Xr / (Xr + Yr + Zr);
y = Yr / (Xr + Yr + Zr);
}
xyy[0] = (int) (255*x + .5);
xyy[1] = (int) (255*y + .5);
xyy[2] = (int) (255*Y + .5);
}
RGB to XYZ
public void rgb2xyz(int R, int G, int B, int []xyz) {
float rf, gf, bf;
float r, g, b, X, Y, Z;
r = R/255.f; //R 0..1
g = G/255.f; //G 0..1
b = B/255.f; //B 0..1
if (r <= 0.04045)
r = r/12;
else
r = (float) Math.pow((r+0.055)/1.055,2.4);
if (g <= 0.04045)
g = g/12;
else
g = (float) Math.pow((g+0.055)/1.055,2.4);
if (b <= 0.04045)
b = b/12;
else
b = (float) Math.pow((b+0.055)/1.055,2.4);
X = 0.436052025f*r + 0.385081593f*g + 0.143087414f *b;
Y = 0.222491598f*r + 0.71688606f *g + 0.060621486f *b;
Z = 0.013929122f*r + 0.097097002f*g + 0.71418547f *b;
xyz[1] = (int) (255*Y + .5);
xyz[0] = (int) (255*X + .5);
xyz[2] = (int) (255*Z + .5);
}
RGB to LAB
public void rgb2lab(int R, int G, int B, int []lab) {
//http://www.brucelindbloom.com
float r, g, b, X, Y, Z, fx, fy, fz, xr, yr, zr;
float Ls, as, bs;
float eps = 216.f/24389.f;
float k = 24389.f/27.f;
float Xr = 0.964221f; // reference white D50
float Yr = 1.0f;
float Zr = 0.825211f;
// RGB to XYZ
r = R/255.f; //R 0..1
g = G/255.f; //G 0..1
b = B/255.f; //B 0..1
// assuming sRGB (D65)
if (r <= 0.04045)
r = r/12;
else
r = (float) Math.pow((r+0.055)/1.055,2.4);
if (g <= 0.04045)
g = g/12;
else
g = (float) Math.pow((g+0.055)/1.055,2.4);
if (b <= 0.04045)
b = b/12;
else
b = (float) Math.pow((b+0.055)/1.055,2.4);
X = 0.436052025f*r + 0.385081593f*g + 0.143087414f *b;
Y = 0.222491598f*r + 0.71688606f *g + 0.060621486f *b;
Z = 0.013929122f*r + 0.097097002f*g + 0.71418547f *b;
// XYZ to Lab
xr = X/Xr;
yr = Y/Yr;
zr = Z/Zr;
if ( xr > eps )
fx = (float) Math.pow(xr, 1/3.);
else
fx = (float) ((k * xr + 16.) / 116.);
if ( yr > eps )
fy = (float) Math.pow(yr, 1/3.);
else
fy = (float) ((k * yr + 16.) / 116.);
if ( zr > eps )
fz = (float) Math.pow(zr, 1/3.);
else
fz = (float) ((k * zr + 16.) / 116);
Ls = ( 116 * fy ) - 16;
as = 500*(fx-fy);
bs = 200*(fy-fz);
lab[0] = (int) (2.55*Ls + .5);
lab[1] = (int) (as + .5);
lab[2] = (int) (bs + .5);
}
RGB to LUV
public void rgb2luv(int R, int G, int B, int []luv) {
//http://www.brucelindbloom.com
float rf, gf, bf;
float r, g, b, X_, Y_, Z_, X, Y, Z, fx, fy, fz, xr, yr, zr;
float L;
float eps = 216.f/24389.f;
float k = 24389.f/27.f;
float Xr = 0.964221f; // reference white D50
float Yr = 1.0f;
float Zr = 0.825211f;
// RGB to XYZ
r = R/255.f; //R 0..1
g = G/255.f; //G 0..1
b = B/255.f; //B 0..1
// assuming sRGB (D65)
if (r <= 0.04045)
r = r/12;
else
r = (float) Math.pow((r+0.055)/1.055,2.4);
if (g <= 0.04045)
g = g/12;
else
g = (float) Math.pow((g+0.055)/1.055,2.4);
if (b <= 0.04045)
b = b/12;
else
b = (float) Math.pow((b+0.055)/1.055,2.4);
X = 0.436052025f*r + 0.385081593f*g + 0.143087414f *b;
Y = 0.222491598f*r + 0.71688606f *g + 0.060621486f *b;
Z = 0.013929122f*r + 0.097097002f*g + 0.71418547f *b;
// XYZ to Luv
float u, v, u_, v_, ur_, vr_;
u_ = 4*X / (X + 15*Y + 3*Z);
v_ = 9*Y / (X + 15*Y + 3*Z);
ur_ = 4*Xr / (Xr + 15*Yr + 3*Zr);
vr_ = 9*Yr / (Xr + 15*Yr + 3*Zr);
yr = Y/Yr;
if ( yr > eps )
L = (float) (116*Math.pow(yr, 1/3.) - 16);
else
L = k * yr;
u = 13*L*(u_ -ur_);
v = 13*L*(v_ -vr_);
luv[0] = (int) (2.55*L + .5);
luv[1] = (int) (u + .5);
luv[2] = (int) (v + .5);
}