www.pudn.com > Imagepro.rar > ImgProlib.cpp
// ******************************************************************************
// * 图象处理算法库
// * ----------------------------------------------------------------------------
// * C++ File: imgprolib.cpp
// * Header File: imgProlib.h
// *
// * Description:
// * The class of "IMAGE2D" has several operations curently:
// *
// * ReadBMP() -- read the binary BMP 2D image file
// *
// * GaussSmooth() -- smooth the image with Gaussian smoothing filter高斯平滑
// * FlattenHistogram() -- equalize the hisrogram of the input image直方图均衡
// *
// * Edge Detector
// * SobelEdgeDetector() -- Edge detection with Sobel operator Sobel边缘检测
// * CannyEdgeDetector() -- The first drivative of Guaussian edge detector Canny边缘检测
// *
// * Thresholding
// * SimpleThresholding() -- One threshold is given as an input parameter 指定单阈值
// * SimpleDoubleThresholding() -- Up and down thresholds are given as two input parameters 指定双阈值
// * GlobalThresholding() -- Threshold is selected by top P% of the pixels 全局阈值法
// *
// * Thinning() -- Thin the edges after thresholding 细化
// * HoughTransform() -- Hough transform with polar form of the line Hough变换
// * LineFitting() -- Line fitting (Peak fitting) 边线拟合
// *
// * Image tools
// * ClearImage() -- Clear image to be 0 清除Image
// * NegativeImage() -- Make the image to be negative image 图像反相
// * FlipImage() -- Flip the image along the X axis 图像翻转
// * CopyImage() -- Copy one image to another image 图像复制
// * BlendImage() -- Combine two images into one image 图像叠合
// ******************************************************************************
#include "stdafx.h"
#include
#include
#include
#include
#include
#include
#include "imgprolib.h"
//-------------------------------------------------------------
// ReadBMP(): reads the binary BMP 2D image file
//-------------------------------------------------------------
BOOL IMAGE2D::ReadBMP(char *filename)
{
FILE *fp;
if ((fp = fopen(filename, "r")) == NULL) return FALSE;
fclose(fp);
AUX_RGBImageRec *image;
image = auxDIBImageLoad(filename);
width = image->sizeX;
height = image->sizeY;
gray_level = 255; // maximum gray value
input_data = new unsigned char[width*height];
output_data = new unsigned char[width*height];
int format = 3;
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
{
int offset = i * width * format + j * format;
input_data[i * width + j]
= (image->data[offset] + image->data[offset + 1]
+ image->data[offset + 2])/3;
}
delete image;
return TRUE;
}
//---------------------------------------------------------------
// FlattenHistogram(): equalize the hisrogram of the input image
//---------------------------------------------------------------
void IMAGE2D::FlattenHistogram(unsigned char * input_data)
{
GenHistogram(input_data);
int *flatten_table = new int[gray_level + 1];
int i, sum_freq = 0;
// create the transformation table of flattened histogram
for(i = 0; i <= gray_level; i ++)
{
sum_freq += histogram[i];
double cc = (double)(gray_level * sum_freq) / (double)(width * height);
flatten_table[i] = (int)(cc);
}
//transform the original image to equalized image
for(i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
output_data[i*width + j] = flatten_table[input_data[i*width + j]];
delete flatten_table;
}
//----------------------------------------------------------
//GenHistogram(): generate the hisrogram of the input image
//----------------------------------------------------------
void IMAGE2D::GenHistogram(unsigned char *inData)
{
histogram = new unsigned int[gray_level + 1];
for(int i = 0; i <= gray_level; i ++)
histogram[i] = 0; // initialized to be 0
for(i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
histogram[inData[i * width + j]] ++;
}
//----------------------------------------------------------------------------------
//GetHistogram(): generate the hisrogram of the input image with size width x height
//----------------------------------------------------------------------------------
unsigned int * IMAGE2D::GetHistogram(unsigned int *inData, int in_width, int in_height)
{
unsigned int * hist = new unsigned int[gray_level + 1];
for(int i = 0; i <= gray_level; i ++)
hist[i] = 0; // initialized to be 0
for(i = 0; i < in_height; i ++)
for(int j = 0; j < in_width; j ++)
{
int k = inData[i * in_width + j];
hist[k] ++;
}
return hist;
}
//---------------------------------------------------------------------
//GaussSmooth(): smooth the input image with Gaussian smoothing filter
//---------------------------------------------------------------------
void IMAGE2D::GaussSmooth(double sigma, int mask_width, unsigned char * input_data,
unsigned char * output_data)
{
// preprocess the mask_width（在此处针对高斯模板，其矩阵为3×3的。）
if(mask_width < 0) mask_width = - mask_width; // mask_width must be positive
if(!(mask_width % 2))mask_width ++; // mask_width must be odd when even
// allocate the 1D mask matrix and computate its value
double *mask = new double[mask_width];
Gauss1D(sigma, mask_width, mask);
unsigned char * tmp_data = new unsigned char[width*height]; // storing horizontal smoothing
int half_mask_width = (mask_width - 1) / 2;
int tmp_num = 3 * half_mask_width; //for boundary
unsigned char * tmp = new unsigned char[tmp_num];
#define _CalBoundary(output_data) \
{ \
double sum_pixel = 0.0; \
for(k = 0; k < mask_width; k ++) \
sum_pixel += mask[k] * tmp[k]; \
_DT(output_data, i, j) =(int)(sum_pixel); \
}
//-------------- CONVOLUTION WITH THE HORIZONAL MASK -----------------
// calculate X-direction smoothed values of the main body of the image
for(int i = 0; i < height; i ++)
for(int j = half_mask_width; j < width - half_mask_width; j ++)
{
double sum_pixel = 0.0;
for(int k = 0; k < mask_width; k ++)
sum_pixel += mask[k] * _DT(input_data, i, j + k - half_mask_width);
_DT(tmp_data, i, j) =(int)(sum_pixel);
}
// calculate the x-direction boundary smoothed pixels of the image
for(i = 0; i < height; i ++)
{
// boundary ( 0 -> half_mask_width )
//----- calculate the tempotary array
for(int k = 0; k < half_mask_width; k ++)
tmp[k] = _DT(input_data, i, width - half_mask_width + k);
for(k = half_mask_width; k < tmp_num; k ++)
tmp[k] = _DT(input_data, i, k - half_mask_width);
//----- claculate the pixel on the left boundary
for(int j = 0; j < half_mask_width; j ++)
_CalBoundary(tmp_data);
// boundary ( (width - half_mask_width) -> (width - 1) )
//----- calculate the tempotary array
int ltmp = 2 * half_mask_width;
for(k = 0; k < ltmp; k ++)
tmp[k] = _DT(input_data, i, width - ltmp);
for(k = ltmp; k < tmp_num; k ++)
tmp[k] = _DT(input_data, i, k - ltmp);
//----- claculate the pixel on the right boundary
for(j = width - half_mask_width; j < width; j ++)
_CalBoundary(tmp_data);
}
//--------------------- END OF CONVOLUTION WITH THE HORIZONAL MASK -------
//-------------- CONVOLUTION WITH THE VERTICAL MASK -----------------
// calculate Y-direction smoothed values of the main body of the image
for(int j = 0; j < width; j ++)
for(int i = half_mask_width; i < height - half_mask_width; i ++)
{
double sum_pixel = 0.0;
for(int k = 0; k < mask_width; k ++)
sum_pixel += mask[k] * _DT(tmp_data, i + k - half_mask_width, j);
_DT(output_data, i, j) =(int)(sum_pixel);
}
// calculate the Y-direction boundary smoothed pixels of the image
for(j = 0; j < width; j ++)
{
// boundary ( Y: 0 -> half_mask_width )
//----- calculate the tempotary array
for(int k = 0; k < half_mask_width; k ++)
tmp[k] = _DT(input_data, height - half_mask_width + k, j);
for(k = half_mask_width; k < tmp_num; k ++)
tmp[k] = _DT(tmp_data, k - half_mask_width, j);
//----- claculate the pixel on the bottom boundary
for(int i = 0; i < half_mask_width; i ++)
_CalBoundary(output_data);
// boundary ( Y: (height - half_mask_width) -> (height - 1) )
//----- calculate the tempotary array
int ltmp = 2 * half_mask_width;
for(k = 0; k < ltmp; k ++)
tmp[k] = _DT(tmp_data, height - ltmp, j);
for(k = ltmp; k < tmp_num; k ++)
tmp[k] = _DT(tmp_data, k - ltmp, j);
//----- claculate the pixel on the top boundary
for(i = height - half_mask_width; i < height; i ++)
_CalBoundary(output_data);
}
//--------------------- END OF CONVOLUTION WITH THE VERTICAL MASK -------
delete tmp_data;
}
//-----------------------------------------------------------
//Gauss1D(): computate the 1D Gauss mask matrix(计算一维高斯模板矩阵)
//-----------------------------------------------------------
void IMAGE2D::Gauss1D(double sigma, int mask_width, double *mask_matrix)
{
double sum_mask = 0.0;
int half_width = (int)(mask_width - 1) / 2;
for(int i = -half_width; i <= half_width; i ++)
{
double gx = exp(-(double)(i * i) / (2.0 * sigma * sigma));
mask_matrix[i + half_width] = gx;
sum_mask += gx;
}
// Normalize the mask matrix
for(i = -half_width; i <= half_width; i ++)
mask_matrix[i + half_width] /= sum_mask;
}
//----------------------------------------------------------------------
// simpleThresholding(): One threshold is given as an input parameter.
//----------------------------------------------------------------------
void IMAGE2D::simpleThresholding(int threshold, unsigned char *Data)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
{
if(_DT(Data, i, j) >= threshold)
{
_DT(Data, i, j) = 255;
}
else
_DT(Data, i, j) = 0;
}
}
//--------------------------------------------------------------------------------------
// simpleDoubleThresholding(): Up and down thresholds are given as two input parameters.
// threshold1: down threshold, threshold2: up threshold
//--------------------------------------------------------------------------------------
void IMAGE2D::simpleDoubleThresholding(unsigned char *Data, int threshold1, int threshold2)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
{
if(_DT(Data, i, j) >= threshold1 && _DT(Data, i, j) <= threshold2)
_DT(Data, i, j) = 255;
else
_DT(Data, i, j) = 0;
}
}
//---------------------------------------------------------------------------
// GlobalThresholding(): Threshold is selected according to top Percent%
// of the pixels.
//---------------------------------------------------------------------------
int IMAGE2D::GlobalThresholding(float percent, unsigned char *Data)
{
int threshold, sum = 0, num = 0;
GenHistogram(Data);
int topPixels = (int)(percent * 0.01 * width * height);
for(int k = gray_level; k >= 0; k --)
{
sum += histogram[k];
if(sum >= topPixels)
{
threshold = k;
break;
}
}
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
{
if(_DT(Data, i, j) >= threshold)
{
_DT(Data, i, j) = 255;
num ++;
}
else
_DT(Data, i, j) = 0;
}
return num; // Returned value
}
//// Begin: Thinning() ==========================================
//---------------------------------------------------------------
// Thinning():
// Hilditch algrithm is the most common way for thinning
// Use 8 neighbor link
//---------------------------------------------------------------
void IMAGE2D::Thinning(unsigned char *inData)
{
int i,j, k;
int s, nrn, xx, yy, x;
unsigned char gray;
unsigned char *tmpData = new unsigned char[height * width];
// Covert the original image data into the value 0 or 1:
// 0--background pixel, 1-- edge pixel
for(i = 0; i < height; i ++)
for(j = 0; j< width; j ++)
{
gray = _DT(inData, i, j);
if(gray == LinePixel) _DT(tmpData, i, j) = 1;
if(gray == BackgdPixel) _DT(tmpData, i, j) = 0;
_DT(inData, i, j) = _DT(tmpData, i, j);
}
do{
s=0;
// if the pixel value is 3 (this pixel needs to be cleared), set the value to be 0
for(i = 0; i < height; i ++)
{
for(j = 0; j < width; j ++)
{
if(_DT(inData, i, j) == 3) _DT(inData, i, j) = 0;
_DT(tmpData, i, j) = _DT(inData, i, j);
}
}
CopyImage(inData, tmpData);
// begin to scan the image
for(i = 1; i < height - 1; i ++)
{
yy = i;
for(j = 1; j < width - 1; j ++)
{
xx = j;
// if pixel isn't 0, don't process it
if(_DT(inData, i, j) != 1) continue;
// if all pixels in 4 neightbor aren't 0, don't process it
if((_DT(tmpData, i-1, j) * _DT(tmpData, i+1, j)
* _DT(tmpData, i, j-1) * _DT(tmpData, i, j+1)) !=0)
continue;
// if the value of arround pixels nrn <=1, set the pixel value 2 (this pixel can't be cleared)
nrn = _DT(tmpData, yy-1, xx-1) + _DT(tmpData, yy-1, xx)
+ _DT(tmpData, yy-1, xx+1) + _DT(tmpData, yy, xx-1)
+ _DT(tmpData, yy, xx+1) + _DT(tmpData, yy+1, xx-1)
+ _DT(tmpData, yy+1, xx) + _DT(tmpData, yy+1, xx+1);
if(nrn <= 1) {
_DT(inData, i, j) = 2; continue;
}
CalCrossNumber(inData, xx, yy, &x, &k);
// if cross number x isn't 1, set the pixel value to be 2
if( x != 1) {
_DT(inData, i, j) = 2; continue;
}
if(_DT(inData, i-1, j) == 3)
{
_DT(inData, i-1, j) = 0;
CalCrossNumber(inData, xx, yy, &x, &k);
// if the cross number isn't 1, set the up neightbor pixel to be 3
if(x != 1)
{
_DT(inData, i-1, j) = 3; continue;
}
_DT(inData, i-1, j) = 3;
}
if(_DT(inData, i, j-1) != 3)
{
_DT(inData, i, j) = 3; s=1; continue;
}
_DT(inData, i, j-1) = 0;
CalCrossNumber(inData, xx, yy, &x, &k);
if(x==1)
{
_DT(inData, i, j-1) = 3;
_DT(inData, i, j) = 3;
s=1;
}
else
_DT(inData, i, j-1) = 3;
}
}
// if have no other pixel needs to be processed, end scanning
if( s==0) break;
}while(1);
for(i=0; i= height - half_sizeY
|| j < half_sizeX || j >= width - half_sizeX)
_DT(outData, i, j) = 0;
else
{
ConvolveOnePixel(sizeX, sizeY, Sx, i, j, &pixelX, inData);
ConvolveOnePixel(sizeX, sizeY, Sy, i, j, &pixelY, inData);
_DT(tmpdata, i, j) = _INT_SQRT(pixelX, pixelY);
if( _DT(tmpdata, i, j) > max)
max = _DT(tmpdata, i, j);
}
}
// ---- scaling
float scale = 1.0f;
if(max > 255) scale = 255.f / (float)max;
else scale = (float)max / 255.f;
for(i = 0; i < height; i ++)
for( int j = 0; j < width; j ++)
_DT(outData, i, j)
= (unsigned char)(_DT(tmpdata, i, j) * scale);
delete tmpdata;
}
//// Begin: CannyEdgeDetector() ==========================================
//--------------------------------------------------------------------
// CannyEdgeDetector(): the first drivative of Guaussian Edge Detector
//
// | -1 1 |
// Sx = | -1 1 |
//
// | 1 1 |
// Sy = | -1 -1 |
//--------------------------------------------------------------------
void IMAGE2D::CannyEdgeDetector(int mask_width, double sigma,
unsigned char *inData, unsigned char *outData)
{
static int Sx[] = {-1, 1, -1, 1};
static int Sy[] = { 1, 1, -1, -1};
static int sizeX = 2, sizeY = 2;
GaussSmooth(sigma, mask_width, inData, outData);
CopyImage(outData, inData);
int *tmpdata = new int[width*height];
int max = 0;
// calculate values of the main body of the image
int pixelX, pixelY;
theta = new int[width * height];
for(int i = 0; i < height - 1; i ++)
for(int j = 0; j < width - 1; j ++)
{
ConvolveOnePixel(sizeX, sizeY, Sx, i, j, &pixelX, inData);
ConvolveOnePixel(sizeX, sizeY, Sy, i, j, &pixelY, inData);
// Compute the magnitude
_DT(outData, i, j) = _INT_SQRT(pixelX, pixelY);
_DT(tmpdata, i, j) = _INT_SQRT(pixelX, pixelY);
if( _DT(tmpdata, i, j) > max)
max = _DT(tmpdata, i, j);
// Compute the orientation
if(pixelX == 0) {
if( pixelY == 0)
_DT(theta, i, j) = 0;
else if(pixelY > 0)
_DT(theta, i, j) = (int)(ZOOMIN * PI / 2.0f);
else
_DT(theta, i, j) = (int)(- ZOOMIN * PI / 2.0f);
}
else
_DT(theta, i, j)
= (int)(ZOOMIN * atanf((float)pixelY/(float)pixelX));
}
// ---- scaling
float scale = 1.0f;
if(max > 255) scale = 255 / (float)max;
for(i = 0; i < height; i ++)
for( int j = 0; j < width; j ++)
_DT(outData, i, j)
= (unsigned char)(_DT(tmpdata, i, j) * scale);
NonmaxSuppression(outData);
int half_size = (mask_width - 1) / 2;
for(i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
// clear to be 0 on the boundary
if(i < half_size || i >= height - half_size
|| j < half_size || j >= width - half_size)
_DT(outData, i, j) = 0;
else continue;
delete tmpdata;
}
//--------------------------------------------------------------------------
// NonmaxSuppression(): Nonmaxima suppresion thins the ridges of gradient
// (Canny algorithms) magnitude by suppressing all values along the line
// of the gradient that are not peak values of a ridge.
//--------------------------------------------------------------------------
void IMAGE2D::NonmaxSuppression(unsigned char * input_data)
{
int i1, i2, j1, j2;
for(int i = 1; i < height - 1; i ++)
for(int j = 1; j < width - 1; j ++)
{
FindOrientNeighbors(i, j, &i1, &j1, &i2, &j2);
if(_DT(input_data, i, j) <= _DT(input_data, i1, j1)
&& _DT(input_data, i, j) <= _DT(input_data, i2, j2))
_DT(input_data, i, j) = 0;
}
delete theta;
}
//--------------------------------------------------------------------------
// FindOrientNeighbors(): Find two neighbors of center pixel along its
// (Canny algorithms) orientation within the 3x3 neighborhood.
//--------------------------------------------------------------------------
void IMAGE2D::FindOrientNeighbors(int i0, int j0, int *i1, int *j1, int *i2, int *j2)
{
float PI_8 = (float)(ZOOMIN * PI) / 8.0f;
float PI_4 = (float)(ZOOMIN * PI) / 4.0f;
float PI_2 = (float)(ZOOMIN * PI) / 2.0f;
float angle = (float)(_DT(theta, i0, j0));
if(angle >= - PI_8 && angle < PI_8
|| angle >= PI - PI_8 && angle < PI + PI_8) {
*i1 = i0 + 1; *j1 = j0;
*i2 = i0 - 1; *j2 = j0;
}
else if(angle >= PI_8 && angle < PI_4 + PI_8
|| angle >= - PI + PI_8 && angle < - PI_2 - PI_8) {
*i1 = i0 + 1; *j1 = j0 + 1;
*i2 = i0 - 1; *j2 = j0 - 1;
}
else if(angle >= PI_2 - PI_8 && angle < PI_2 + PI_8
|| angle >= - PI_2 - PI_8 && angle < - PI_2 + PI_8) {
*i1 = i0; *j1 = j0 + 1;
*i2 = i0; *j2 = j0 - 1;
}
else {
*i1 = i0 - 1; *j1 = j0 + 1;
*i2 = i0 + 1; *j2 = j0 - 1;
}
}
//// End of CannyEdgeDetector() ====================================
//---------------------------------------------------------------
// ConvolveOnePixel(): Convolution for one pixel: ImgData[i,j]
//---------------------------------------------------------------
void IMAGE2D::ConvolveOnePixel(int sizeX, int sizeY, int* mask, int i, int j,
int* pixel_value, unsigned char* inData)
{
int half_sizeX = (sizeX - 1) / 2;
int half_sizeY = (sizeY - 1) / 2;
// calculate values of the main body of the image
int sum = 0;
for(int n = 0; n < sizeY; n ++)
for(int m = 0; m < sizeX; m ++)
sum += mask[n * sizeX + m]
* _DT(inData, i - half_sizeY + n, j - half_sizeX + m);
*pixel_value = sum;
}
//// Begin: HoughTransform() ===========================================
//----------------------------------------------------------------------
// HoughTransform():
// line representation: p = x * cos(theta) + y * sin(theta)
//----------------------------------------------------------------------
void IMAGE2D::HoughTransform(unsigned char *inData, float bins)
{
// Thresholding first (Two values)
//---- For displaying accumulator (rtData) ----------------------------
int roh_num = (int)(height / bins); // the number of divided roh
int theta_num =(int)( width / bins); // the number of divided theta
unsigned char *rtData = new unsigned char[height * width];
ClearImage(rtData);
// set the diagonal of the original image to be the maxium of p
rohmax =2 * (float)sqrt((float) width*width + (float) height*height);
Hough_width = Hough_height = (int)(rohmax / bins + 0.5f);
HoughData = new unsigned int[Hough_height * Hough_width];
// Clear to be 0
for(int i = 0; i < Hough_height; i ++)
for( int j = 0; j < Hough_width; j ++)
HoughData[i * Hough_width + j] = 0;
float theta, roh;
int pre_r = 0, pre_roh = 0, rohi;
for(int y = 0; y < height; y ++)
for(int x = 0; x < width; x ++)
{
// if the pixel value is a background pixel, continue
if(_DT(inData, y, x) == BackgdPixel) continue;
//x*cos(theta)+y*sin(theta)-p=0
//theta=[-PI, PI), p=(-rohmax, rohmax)
// ---- for HoughData --------
for(int t = 0; t < Hough_width; t ++)
{
//theta=[-PI,PI)
theta = (float)((2 * (float)t / Hough_width - 1) * PI);
roh =(float)(x*cos(theta) + y*sin(theta));
rohi = (int)(roh*Hough_height/rohmax + Hough_height/2 + 0.5f);
HoughData[rohi * Hough_width + t] ++;
if(t > 0)
{
int droh = rohi - pre_roh;
int sym = 1;
if(droh<0) sym = -1;
if(abs(droh) > 1)
{
for(int n = 1; n < abs(droh) ; n ++)
if(n <= abs(droh/2))
HoughData[(pre_roh + n*sym) * Hough_width + t - 1] ++;
else
HoughData[(pre_roh + n*sym) * Hough_width + t] ++;
}
}
pre_roh = rohi;
}
//----- for rtData -----
for(t = 0; t < theta_num; t ++)
{
//theta=[-PI,PI)
theta = (float)((2 * (float)t / theta_num - 1) * PI);
roh =(float)(x*cos(theta) + y*sin(theta));
int r = (int)((roh + rohmax/2)*(float)roh_num/rohmax + 0.5f);
_DT(rtData, r, t) ++;
if(t > 0)
{
int dr = r - pre_r;
int sym = 1;
if(dr<0) sym = -1;
if(abs(dr) > 1)
{
for(int n = 1; n < abs(dr) ; n ++)
if(n <= abs(dr/2))
_DT(rtData, pre_r + n*sym, t - 1) ++;
else
_DT(rtData, pre_r + n*sym, t) ++;
}
}
pre_r = r;
}
}
CopyImage(rtData, inData);
delete rtData;
}
//// End of HoughTransform() ===========================================
//// Begin: LineFitting() ==============================================
//----------------------------------------------------------------------
// LineFitting():
//----------------------------------------------------------------------
void IMAGE2D::LineFitting(unsigned char *inData, float thresh_percent)
{
int threshold;
int num_peaks = PeakDetection(inData, thresh_percent, &threshold);
int *peaks = new int[2*num_peaks];
int numOfPeaks = GetPeaks(threshold, num_peaks, peaks);
numOfPeaks = PeakFitting(numOfPeaks, peaks); // new peaks
BackProjection(inData, peaks, numOfPeaks);
}
//----------------------------------------------------------------------
// BackProjection():
//----------------------------------------------------------------------
void IMAGE2D::BackProjection(unsigned char *outData, int *peaks, int numOfPeaks)
{
ClearImage(outData);
// --- Calculate the (x,y) along X-axi
for(int x = 0; x < width; x ++)
for( int n = 0; n < numOfPeaks; n ++)
{
// theta=[-PI, PI)
float theta = (float)(2* PI * (float)peaks[n] / Hough_width - PI);
float roh = ((float)peaks[numOfPeaks + n]/Hough_height -0.5f)* rohmax;
int yy = (int)((roh - x * cos(theta)) / sin(theta));
if(yy >= 0 && yy < height)
_DT(outData, yy, x) = LinePixel;
}
// Calculate the (x,y) along Y-axi
for(int y = 0; y < height; y ++)
for( int n = 0; n < numOfPeaks; n ++)
{
// theta=[-PI, PI)
float theta = (float)(2* PI * (float)peaks[n] / Hough_width - PI);
float roh = ((float)peaks[numOfPeaks + n]/Hough_height -0.5f)* rohmax;
int xx = (int)((roh - y * sin(theta)) / cos(theta));
if(xx >= 0 && xx < width)
_DT(outData, y, xx) = LinePixel;
}
}
//----------------------------------------------------------------------
// PeakFitting(): return the number of new peaks
//----------------------------------------------------------------------
int IMAGE2D::PeakFitting(int numOfPeaks, int *peaks)
{
int dx, dy;
dx = dy = 5;
float rohmax = 2 * (float)sqrt((float) width*width + (float) height*height);
int theta_num = (int)(rohmax + 0.5f);
int *newX = new int[numOfPeaks];
int *newY = new int[numOfPeaks];
newX[0] = peaks[0];
newY[0] =peaks[numOfPeaks];
int newNum = 1, flag;
for(int n = 1; n < numOfPeaks; n ++)
{
int curX = peaks[n];
int curY = peaks[numOfPeaks + n];
flag = 0;
for(int m = 0; m < newNum; m ++)
{
// If similar, calculate their average value
if(abs(curX - newX[m]) <= dx
&& abs(curY - newY[m]) <= dy)
{
if(abs(curX - theta_num/2) < dx) // theta = 0
newX[m] = theta_num/2;
else if(abs(curX - theta_num/4) < dx) // theta = -PI/2
newX[m] = theta_num / 4;
else if(abs(curX - theta_num*3/4) < dx) // theta = PI/2
newX[m] = theta_num * 3 / 4;
else
newX[m] =(int)((newX[m] + curX) / 2.0f + 0.5);
newY[m] = (int)((newY[m] + curY) / 2.0f + 0.5);
flag = 1;
}
}
if(flag == 0) // not found the similar peak
{
newX[newNum] = curX;
newY[newNum] = curY;
newNum ++;
}
}
for(int i = 0; i < newNum; i ++)
{
peaks[i] = newX[i];
peaks[i + newNum] = newY[i];
}
return newNum;
}
//----------------------------------------------------------------------
// GetPeaks():
//----------------------------------------------------------------------
int IMAGE2D::GetPeaks(int threshold, int numOfPeaks, int *peaks)
{
int num = 0;
int width1 = (int)((float)Hough_width/4.0f);
int width2 = (int)((float)Hough_width*3.0f/4.0f);
int * peakX = new int[numOfPeaks];
int * peakY = new int[numOfPeaks];
for(int i = 0; i < Hough_height; i ++)
for(int j = width1; j < width2; j ++) // theta=[-PI/2, PI/2)
if(HoughData[i * Hough_width + j] >= (unsigned int)threshold)
{
peakY[num] = i; // roh
peakX[num] = j; // theta
num ++;
}
numOfPeaks = num;
for( int k = 0; k < numOfPeaks; k ++)
{
peaks[numOfPeaks + k] = peakY[k];
peaks[k] = peakX[k];
}
return numOfPeaks;
}
//----------------------------------------------------------------------
// PeakDetection():
//----------------------------------------------------------------------
int IMAGE2D::PeakDetection(unsigned char * inData, float percent, int *threshold)
{
//Peak detection (Thresholding algorithm similar to global thresholding algoroghtm)
unsigned int * hist = new unsigned int[gray_level + 1];
hist = GetHistogram(HoughData, Hough_width, Hough_height);
int sum = 0;
int topPixels = (int)(percent * 0.01 * ( Hough_width * Hough_height - hist[0]));
for(int k = gray_level; k > 0; k --) // except k = 0
{
sum += hist[k];
if(sum >= topPixels)
{
*threshold = k;
break;
}
}
int numOfPeaks = sum;
return numOfPeaks;
}
//// End of LineFitting() ==============================================
//// Begin: Image tools ===================================
//---------------------------------------------------------
// NegativeImage(): Make the image to be negative image
//---------------------------------------------------------
void IMAGE2D::NegativeImage(unsigned char *data)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
_DT(data, i, j) = gray_level - _DT(data, i, j);
}
//---------------------------------------------
// FlipImage(): Filp the image along the X axis
//---------------------------------------------
void IMAGE2D::FlipImage(unsigned char *data)
{
tmp_data = new unsigned char[width * height];
CopyImage(data, tmp_data);
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
_DT(data, i, j) = _DT(tmp_data, height - i, j);
delete tmp_data;
}
//---------------------------------------------------------
// CopyImage(): Copy one image to another image
//---------------------------------------------------------
void IMAGE2D::CopyImage(unsigned char *srcData, unsigned char *data)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
data[i * width + j] = srcData[i * width + j];
}
//---------------------------------------------------------
// ClearImage(): Clear image to be 0
//---------------------------------------------------------
void IMAGE2D::ClearImage(unsigned char *data)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
data[i*width+j] = 0;
}
//---------------------------------------------------------
// BlendImage(): Combine two images into one image
//---------------------------------------------------------
void IMAGE2D::BlendImage(unsigned char *srcData, unsigned char *data)
{
for(int i = 0; i < height; i ++)
for(int j = 0; j < width; j ++)
data[i*width+j] |= srcData[i*width+j];
}
// ************************************************************************
//
// 图像形态学变换函数
//
// ErosionGrayImage() - 图像腐蚀
// DilationGrayImage() - 图像膨胀
// OpenGrayImage() - 图像开运算
// CloseGrayImage() - 图像闭运算
//
// ************************************************************************
/*************************************************************************
*
* 函数名称：
* ErosionGrayImage()
*
* 参数:
* const unsigned char *srcData - 指向源图像指针
* unsigned char *dstData - 指向目标图像的指针
* int element[2][2] - 结构元数组
*
* 返回值:
* BOOL - 腐蚀成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用于对图像进行腐蚀运算。
* 目标图像为灰度图像。
************************************************************************/
BOOL IMAGE2D::ErosionGrayImage(unsigned char *srcData, unsigned char *dstData, int element[2][2])
{
int i,j;
int m=height,n=width;
unsigned char gray;
for(i=0;it)
j--;
x[i]=x[j];
while (i graylevel2)
{
w = 2 * (graylevel1 - graylevel2);
}
else
{
w = 2 * (graylevel2 - graylevel1);
}
if(h_fre1 == h && h_fre2 == h)
{
m = 1;
}
else
{
m = 2;
}
}
void IMAGE2D::GradientAmplitudeHistogram(unsigned char *WndData)//产生梯度幅度直方图
{
int i,j;
//计算图像上各灰度级对应的灰度梯度幅度值
for(i = 0; i < ScoutWndWidth; i++)
{
for(j = 0; j < ScoutWndHeight; j++)
{
gray[i * ScoutWndWidth + j] = WndData[i * ScoutWndWidth +j];
GradsAmp[gray[i * width + j]] =
sqrt((gray[i * ScoutWndWidth + j] - gray[(i + 1) * ScoutWndWidth+ j])
* (gray[i * ScoutWndWidth + j] - gray[(i + 1) * ScoutWndWidth + j])
+ (gray[i * ScoutWndWidth + j] - gray[i * ScoutWndWidth + j + 1])
* (gray[i *ScoutWndWidth + j] - gray[i * ScoutWndWidth + j + 1]));
}
}
//产生灰度梯度幅度直方图（统计各灰度级值的梯度幅度（大小）的分布情况）
//计算背景的灰度梯度幅度均值和均方差
double sum_GradsAmp1 = sum_element(GradsAmp,gray_level1);
for(i = 0;i < gray_level1; i++)
{
pr1[i] = (GradsAmp[i])/(sum_GradsAmp1);
part1[i] = GrayLevel[i] * pr1[i];
}
m_GA1 = sum_element(part1,gray_level1);
for(i = 0; i < gray_level1; i++)
{
part_MSE1[i] = (part1[i] - m_GA1) * (part1[i] - m_GA1);
}
double sum_MSEGA1 = sum_element(part_MSE1,gray_level1);
MSE_GA1 = sqrt(sum_MSEGA1/(gray_level1 - 1));
//计算目标的灰度梯度幅度均值和均方差
int n = graylevel - gray_level2;
double sum_GradsAmp2 = sum_element(GradsAmp,n);
for(i = gray_level2;i < graylevel; i++)
{
pr1[i] = (GradsAmp[i])/(sum_GradsAmp2);
part1[i] = GrayLevel[i] * pr1[i];
}
m_GA2 = sum_element(part1,n);
for(i = gray_level2; i < graylevel; i++)
{
part_MSE1[i] = (part1[i] - m_GA2) * (part1[i] - m_GA2);
}
double sum_MSEGA2 = sum_element(part_MSE1,n);
MSE_GA2 = sqrt(sum_MSEGA2/(n - 1));
//generate the gray gradient amplitude histogram of the input image
for(i = 0; i < graylevel; i ++)
{
x[i] = -(i-m_GA1)*(i-m_GA1)/((MSE_GA1*MSE_GA1)*2);
Math_GA1 = p1/(MSE_GA1*sqrt(2*PI));
GrayGA_fre1[i] = Math_GA1*(exp(x[i]));
y[i] = -(i-m_GA2)*(i-m_GA2)/((MSE_GA2*MSE_GA2)*2);
Math_GA2 = p2/(MSE_GA2*sqrt(2*PI));
GrayGA_fre2[i] = Math_GA2*(exp(x[i]));
//产生高斯拟合灰度梯度幅度直方图
GradsAmp[i] = GrayGA_fre1[i] + GrayGA_fre2[i];
}
//梯度幅度直方图的峰高和峰宽
h1_fre1 = max_element(GrayGA_fre1,gray_level1);
h1_fre2 = max_element(GrayGA_fre2,n);
h1 = max_element(GradsAmp,graylevel);
graylevel11 = graylevel;
min_element(GradsAmp,gray_level);
graylevel12 = graylevel;
if(graylevel11 > graylevel12)
{
w = 2 * (graylevel11 - graylevel12);
}
else
{
w = 2 * (graylevel12 - graylevel11);
}
if(h_fre1 == h && h_fre2 == h)
{
m = 1;
}
else
{
m = 2;
}
}
void IMAGE2D::GradientDirectionHistogram(unsigned char *WndData)
{
int i,j;
for(i = 0; i < ScoutWndWidth; i++)
{
for(j = 0; j < ScoutWndHeight; j++)
{
gray[i * ScoutWndWidth + j] = WndData[i * ScoutWndWidth +j];
GradsDirect[gray[i * ScoutWndWidth + j]] =
(gray[i * ScoutWndWidth + j] - gray[i * ScoutWndWidth + j + 1])
/(gray[i * ScoutWndWidth + j] - gray[(i + 1) * ScoutWndWidth + j]);
}
}
//产生灰度梯度方向直方图（统计各灰度级值的梯度方向的分布情况）
//计算背景的灰度梯度方向均值和均方差
double sum_GradsDirect1 = sum_element(GradsDirect,gray_level1);
for(i = 0;i < gray_level1; i++)
{
pr2[i] = (GradsDirect[i])/(sum_GradsDirect1);
part2[i] = GrayLevel[i] * pr2[i];
}
m_GD1 = sum_element(part2,gray_level1);
for(i = 0; i < gray_level1; i++)
{
part_MSE2[i] = (part2[i] - m_GD1) * (part2[i] - m_GD1);
}
double sum_MSEGD1 = sum_element(part_MSE2,gray_level1);
MSE_GD1 = sqrt(sum_MSEGD1/(gray_level1 - 1));
//计算目标的灰度梯度方向均值和均方差
int n = graylevel - gray_level2;
double sum_GradsDirect2 = sum_element(GradsDirect,n);
for(i = gray_level2;i < graylevel; i++)
{
pr2[i] = (GradsDirect[i])/(sum_GradsDirect2);
part2[i] = GrayLevel[i] * pr2[i];
}
m_GD2 = sum_element(part2,n);
for(i = gray_level2; i < graylevel; i++)
{
part_MSE2[i] = (part2[i] - m_GD2) * (part2[i] - m_GD2);
}
double sum_MSEGD2 = sum_element(part_MSE2,n);
MSE_GD2 = sqrt(sum_MSEGD2/(n - 1));
MSE_GD = MSE_GD1 + MSE_GD2;
// a1' = 1/MSE_GD;
//generate the gray gradient amplitude histogram of the input image
for(i = 0; i < graylevel; i ++)
{
x[i] = -(i-m_GD1)*(i-m_GD1)/((MSE_GD1*MSE_GD1)*2);
Math_GD1 = p1/(MSE_GD1*sqrt(2*PI));
GrayGD_fre1[i] = Math_GD1*(exp(x[i]));
y[i] = -(i-m_GD2)*(i-m_GD2)/((MSE_GD2*MSE_GD2)*2);
Math_GD2 = p2/(MSE_GD2*sqrt(2*PI));
GrayGD_fre2[i] = Math_GD2*(exp(x[i]));
//产生高斯拟合灰度梯度方向直方图
GradsDirect[i] = GrayGD_fre1[i] + GrayGD_fre2[i];
}
//梯度方向直方图的峰高和峰宽
h2_fre1 = max_element(GrayGD_fre1,gray_level1);
h2_fre2 = max_element(GrayGD_fre2,n);
h2 = max_element(GradsDirect,graylevel);
graylevel21 = graylevel;
min_element(GradsDirect,graylevel);
graylevel22 = graylevel;
if(graylevel21 > graylevel22)
{
w = 2 * (graylevel21 - graylevel22);
}
else
{
w = 2 * (graylevel22 - graylevel21);
}
if(h_fre1 == h && h_fre2 == h)
{
m = 1;
}
else
{
m = 2;
// if(w == 16)
// {
// a2' =1;
//}
// else if(w > 10 &&w !=16)
// {
// a2' = 0;
// }
//else
//{
// a2' = 0.5;
//}
}
}
void IMAGE2D::BasicDollarScout(unsigned char *inData,int BasDolSize)
{
int j = BasDolSize;
//定义基元搜索窗口的大小，并指定其初始位置
if(m_Parms.fMakeUserHitOKToScout == TRUE)
m_Parms.iWindowHeight = BasDolSize + 5;//设置窗口的高度
m_Parms.iWindowWidth = 30;//inpixels，设置窗口的宽度
m_ScoStatus.x = 0;
m_ScoStatus.y = 0;//窗口的初始位置
m_Parms.diScout = eLEFTtoRIGHT;
int mNumberMax = (int)((width*height)/(m_Parms.iWindowWidth*m_Parms.iWindowHeight));//搜索窗口的最大值
//搜索窗口在影像上滑动，从第一个窗口开始道路基元的搜索
for(int i = 1; i <= mNumberMax; i++)
{
m_ScoStatus.x += 30;
m_ScoStatus.y += (BasDolSize + 5);
if(m_ScoStatus.fBeginScout = TRUE)
{
//判断梯度方向直方图是否存在对称双峰
for(; j >= BasDolSize - 10;j--)
{
for(int x = m_ScoStatus.x; x <= m_ScoStatus.x +30; x++)
for(int y = m_ScoStatus.y; y <= m_ScoStatus.y+j; y++)
{
IMAGE2D::GradientDirectionHistogram(DolData);
if(m == 1)
{
return;
}
else if(m == 2)
{
IMAGE2D::GaussFittingHistogram(DolData);
if(m == 2)
{
return;
}
else if(m == 1)
{
IMAGE2D::GradientAmplitudeHistogram(DolData);
if(m == 2)
{
return;
}
else if(m ==1 && h != 0)
{
m_ScoStatus.fScoutBasicDollar = TRUE;
}
}
}
}
}
}
}
}
void IMAGE2D::RoadEdgeAffirm(unsigned char *SlipWinData)
{
}
BOOL IMAGE2D::Write(CFile *pFile)
{
BITMAPFILEHEADER bmfh;
// 设置文件头中文件类型为"BM"
bmfh.bfType = 0x4d42;
// 计算信息头和调色板的大小尺寸
int nSizeHdr = sizeof(BITMAPINFOHEADER) + sizeof(RGBQUAD) * m_nColorTableEntries;
// 指定文件的大小
bmfh.bfSize = sizeof(BITMAPFILEHEADER) + nSizeHdr + m_dwSizeImage;
//设置保留字
bmfh.bfReserved1 = bmfh.bfReserved2 = 0;
//从文件头到实际的位图数据的偏移字节数（除了位图数据以外的前三部分之和）
bmfh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) +
sizeof(RGBQUAD) * m_nColorTableEntries;
// 进行写操作
try
{
//文件头的写入
pFile->Write((LPVOID) &bmfh, sizeof(BITMAPFILEHEADER));
//信息头和调色板的写入
pFile->Write((LPVOID) m_lpBMIH, nSizeHdr);
//实际位图数据的写入
pFile->Write((LPVOID) m_lpImage, m_dwSizeImage);
}
// 错误处理
catch(CException* pe)
{
pe->Delete();
AfxMessageBox("write error");
return FALSE;
}
// 返回
return TRUE;
}