利用梯度方向在图像上提取直线的方法原始文献：J.Brian Burns, Allen R.Hanson,Edward M.Riseman, Extracti
利用梯度方向在图像上提取直线的方法
原始文献：
J.Brian Burns, Allen R.Hanson,Edward M.Riseman, Extracting Straight Lines,IEEE Transactions on Pattern Analysis and Machine Intelligence,(Volume:PAMI-8 , Issue: 4 ),1986
文章下载地址：http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4767808&tag=1，如果无法下载IEEE的文章就去我的资源里下载：http://download.csdn.net/detail/ihadl/6446169
作者主页：http://www.ai.sri.com/people/burns/，这里虽然列了这篇文章，但是没有找到源代码，只有一些作者的软件（http://www.ai.sri.com/software_list/），看着都好高端的样子，以后有时间可以下下来看看效果。
在pudn上下到一个该文章的源代码（http://www.pudn.com/downloads204/sourcecode/graph/detail960709.html），不知道是原作者实现的还是后人实现的，总之最后历尽千辛万苦还是可以调试成功的，也是可用的，在此感谢代码实现者！
然后就对该算法进行了研究，对程序进行了调试，前后花费近一个月时间，因此有必要总结一下。
算法的核心步骤：Introduction中介绍了一些提取线的通用方法和存在的一些问题，其中针对的主要问题是大多数算法利用了像素梯度的大小，而没有利用梯度的方向，因此文章主要利用梯度的方向来提取直线。主要四个步骤：
一、Group pixels into line-support regions based on similarity of gradient orientation
根据梯度的方向构造LSR.主要步骤包括:
(1)选择梯度算子计算x和y方向的梯度,然后计算梯度的方向角度.
(2)利用固定间隔的区间对梯度角度结果进行分割.由于区域增长分割方法产生的偶然错误都会对结果产生巨大影响.因此使用一种固定分区的分割方法.首先将角度360度的角度范围分割为4个区间或者8个区间,然后将具体的角度归并入分好的角度区间,最后再利用邻域搜索的算法对角度区间的结果进行聚合.
(3)利用重叠间隔的区间对梯度角度结果进行分割.由于利用固定的区间会产生问题(比如不同的两条直线由于相邻且角度相近被聚合为一条;一条直线很可能被分到两个不同的LSR当分割区间恰巧将它分开).因此扩展到重叠区间的方法,先以0度为起点,45度的间隔进行一次归并,再以22.5度为起点,45度的间隔进行第二次归并,然后利用获取最长线的原则将两次的结果合并.具体步骤是:1)先获取每个LSR的长度.2)如果某个像素包含在两个不同的LSR中,它就投票给长度长的那个LSR.3)每个LSR都会得到内部像素的投票.4)取投票数在50%以上的作为LSR.大部分的LSR要么获得非常多的票,要么非常少.
二.Approximate the intensity surface by a planar surface
获得LSR以后需要确定准确的直线的位置。第一个先验知识是直线的方向垂直于LSR中各点的梯度方向。第二就是确定沿着这个直线方向的直线的具体位置。利用两平面横切的方法获取具体位置。第一个平面是梯度大小利用最小二乘拟合得到的平面，第二个平面是图像原始灰度值的平均得到的水平平面，这两个平面相交就获得了最终的直线位置。
三.Extract attributes from the line-support region and planar fit.
提取到直线以后可以计算一些线特征，文中提供了一些线特征的计算公式。
四.Filter lines on the attributes to isolate various image eventssuch as long straight lines of any contrast.
根据线特征的一些先验知识剔除不想要的线，留下目标直线。
源代码
具体代码由三个文件组成，首先是bitmap.h和bitmap.c，看代码说明好像是发现了windows.h中的bitmap类有问题，所以就自己重写的bmp的数据结构和读写函数
//
// bitmap.h
//
// header file for MS bitmap format
//
//
#ifndef BITMAP_H
#define BITMAP_H
unsigned char* read_bmp(const char* fname, int* width, int* height, int* step);
void write_bmp(const char *iname, int width, int height, unsigned char *data);
#endif
//
// bitmap.c
//
// handle MS bitmap I/O. For portability, we don't use the data structure defined in Windows.h
// However, there is some strange thing, the side of our structure is different from what it
// should though we define it in the same way as MS did. So, there is a hack, we use the hardcoded
// constanr, 14, instead of the sizeof to calculate the size of the structure.
// You are not supposed to worry about this part. However, I will appreciate if you find out the
// reason and let me know. Thanks.
//
#include "bitmap.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define BMP_BI_RGB 0L
typedef unsigned short BMP_WORD;
typedef unsigned int BMP_DWORD;
typedef int BMP_LONG;
typedef struct {
BMP_WORD bfType;
BMP_DWORD bfSize;
BMP_WORD bfReserved1;
BMP_WORD bfReserved2;
BMP_DWORD bfOffBits;
} BMP_BITMAPFILEHEADER;
typedef struct {
BMP_DWORD biSize;
BMP_LONG biWidth;
BMP_LONG biHeight;
BMP_WORD biPlanes;
BMP_WORD biBitCount;
BMP_DWORD biCompression;
BMP_DWORD biSizeImage;
BMP_LONG biXPelsPerMeter;
BMP_LONG biYPelsPerMeter;
BMP_DWORD biClrUsed;
BMP_DWORD biClrImportant;
} BMP_BITMAPINFOHEADER;
BMP_BITMAPFILEHEADER bmfh;
BMP_BITMAPINFOHEADER bmih;
unsigned char* read_bmp(const char *fname, int* width, int* height, int* step) {
FILE* file=fopen(fname, "rb");
if (!file) return NULL;
// I am doing fread( &bmfh, sizeof(BMP_BITMAPFILEHEADER), 1, file ) in a safe way.
fread( &(bmfh.bfType), 2, 1, file);
fread( &(bmfh.bfSize), 4, 1, file);
fread( &(bmfh.bfReserved1), 2, 1, file);
fread( &(bmfh.bfReserved2), 2, 1, file);
fread( &(bmfh.bfOffBits), 4, 1, file);
BMP_DWORD pos = bmfh.bfOffBits;
fread( &bmih, sizeof(BMP_BITMAPINFOHEADER), 1, file );
// error checking
// "BM" actually
if ( bmfh.bfType!= 0x4d42) return NULL;
if ( bmih.biBitCount != 24) return NULL;
fseek( file, pos, SEEK_SET );
*width = bmih.biWidth;
*height = bmih.biHeight;
int padWidth = *width * 3;
int pad = 0;
if ( padWidth % 4!= 0) {
pad = 4- (padWidth % 4);
padWidth += pad;
}
*step = padWidth;
int bytes = *height * padWidth;
unsigned char* data = malloc(bytes);
int foo = fread( data, bytes, 1, file );
if (!foo) {
free(data);
return NULL;
}
fclose( file );
// shuffle bitmap data such that it is (R,G,B) tuples in row-major order
int i, j;
j = 0;
unsigned char temp;
unsigned char* in;
unsigned char* out;
in = data;
out = data;
for ( j = 0; j < *height; ++j ) {
for ( i = 0; i < *width; ++i ) {
out[1] = in[1];
temp = in[2];
out[2] = in[0];
out[0] = temp;
in += 3;
out += 3;
}
in += pad;
}
return data;
}
void write_bmp(const char *iname, int width, int height, unsigned char *data) {
int bytes, pad;
bytes = width * 3;
pad = (bytes%4) ? 4-(bytes%4) : 0;
bytes += pad;
bytes *= height;
bmfh.bfType = 0x4d42; // "BM"
bmfh.bfSize = sizeof(BMP_BITMAPFILEHEADER) + sizeof(BMP_BITMAPINFOHEADER) + bytes;
bmfh.bfReserved1 = 0;
bmfh.bfReserved2 = 0;
bmfh.bfOffBits = /*hack sizeof(BMP_BITMAPFILEHEADER)=14, sizeof doesn't work?*/
14 + sizeof(BMP_BITMAPINFOHEADER);
bmih.biSize = sizeof(BMP_BITMAPINFOHEADER);
bmih.biWidth = width;
bmih.biHeight = height;
bmih.biPlanes = 1;
bmih.biBitCount = 24;
bmih.biCompression = BMP_BI_RGB;
bmih.biSizeImage = 0;
bmih.biXPelsPerMeter = (int)(100/ 2.54* 72);
bmih.biYPelsPerMeter = (int)(100/ 2.54* 72);
bmih.biClrUsed = 0;
bmih.biClrImportant = 0;
FILE *foo=fopen(iname, "wb");
// fwrite(&bmfh, sizeof(BMP_BITMAPFILEHEADER), 1, foo);
fwrite( &(bmfh.bfType), 2, 1, foo);
fwrite( &(bmfh.bfSize), 4, 1, foo);
fwrite( &(bmfh.bfReserved1), 2, 1, foo);
fwrite( &(bmfh.bfReserved2), 2, 1, foo);
fwrite( &(bmfh.bfOffBits), 4, 1, foo);
fwrite(&bmih, sizeof(BMP_BITMAPINFOHEADER), 1, foo);
bytes /= height;
unsigned char* scanline = malloc(bytes);
int i, j;
for (j = 0; j < height; ++j ) {
memcpy( scanline, data + j*3*width, bytes );
for (i = 0; i < width; ++i ) {
unsigned char temp = scanline[i*3];
scanline[i*3] = scanline[i*3+2];
scanline[i*3+2] = temp;
}
fwrite( scanline, bytes, 1, foo);
}
free(scanline);
fclose(foo);
}
接下来是算法的实现函数，burns.h和burns.c,其中使用了 Intel 的IPP图像库（Integrated Performance Primitives），该库的下载地址在后面列出
//
// burns.h
//
#ifndef _BURNS_H
#define _BURNS_H
#include <ipp.h>
/**
* Straight line segment.
*/
typedef struct Line {
/** Coordinates of endpoints. */
double x1, y1, x2, y2;
/** Length of line. */
double length;
/** Number of pixels in the line's support. */
unsigned int num_pixels;
/** Number of pixels who voted for line. */
unsigned int num_votes;
/** Slope in the x-direction. */
double tangent_x;
/** Slope in the y-direction. */
double tangent_y;
/** Midpoint in the x-direction. */
double midpoint_x;
/** Midpoint in the y-direction. */
double midpoint_y;
} Line;
void print_lines(Line* line, int num_lines);
void free_lines(Line* line);
/**
* Find straight line segments in a grayscale image.
* @Article{ Burns1,
* title = "Extracting Straight Lines",
* author = "J. B. Burns and A. R. Hanson and E. M. Riseman",
* journal = "IEEE Transactions on Pattern Analysis and Machine Intelligence",
* pages = "425--455",
* month = jul,
* year = "1986"
* }
*
* @param pixels Array of bytes, the image data.
* @param step Interval in bytes between consecutive rows in the image.
* @param width Width of region of interest.
* @param height Height of region of interest.
* @param roi Region of interest, the size of the image.
* @param num_buckets Increasing the number of buckets decreases the tolerance
* for variations in the gradient orientation along a line. Recommended
* value: 8.
* @param min_gradient Treat gradients below this magnitude (taxi cab distance [0..255])
* as zero. Increasing this value eliminates lines with fuzzy boundaries.
* @param min_length Do not return lines below this length.
*/
void burns_line_extraction(Ipp8u *image, int step, int width, int height,
int num_buckets, int min_gradient, double min_length,
Line** lines, int* num_lines);
void burns_line_extraction_demo(Line** lines, int* num_lines);
#endif /* _BURNS_H */
#ifndef _BURNS_H
#define _BURNS_H
#include <ipp.h>
/**
* Straight line segment.
*/
typedef struct Line {
/** Coordinates of endpoints. */
double x1, y1, x2, y2;
/** Length of line. */
double length;
/** Number of pixels in the line's support. */
unsigned int num_pixels;
/** Number of pixels who voted for line. */
unsigned int num_votes;
/** Slope in the x-direction. */
double tangent_x;
/** Slope in the y-direction. */
double tangent_y;
/** Midpoint in the x-direction. */
double midpoint_x;
/** Midpoint in the y-direction. */
double midpoint_y;
} Line;
void print_lines(Line* line, int num_lines);
void free_lines(Line* line);
/**
* Find straight line segments in a grayscale image.
* @Article{ Burns1,
* title = "Extracting Straight Lines",
* author = "J. B. Burns and A. R. Hanson and E. M. Riseman",
* journal = "IEEE Transactions on Pattern Analysis and Machine Intelligence",
* pages = "425--455",
* month = jul,
* year = "1986"
* }
*
* @param pixels Array of bytes, the image data.
* @param step Interval in bytes between consecutive rows in the image.
* @param width Width of region of interest.
* @param height Height of region of interest.
* @param roi Region of interest, the size of the image.
* @param num_buckets Increasing the number of buckets decreases the tolerance
* for variations in the gradient orientation along a line. Recommended
* value: 8.
* @param min_gradient Treat gradients below this magnitude (taxi cab distance [0..255])
* as zero. Increasing this value eliminates lines with fuzzy boundaries.
* @param min_length Do not return lines below this length.
*/
void burns_line_extraction(Ipp8u *image, int step, int width, int height,
int num_buckets, int min_gradient, double min_length,
Line** lines, int* num_lines);
void burns_line_extraction_demo(Line** lines, int* num_lines);
#endif /* _BURNS_H */
burns.c
//
// burns.c
//
#include "burns.h"
#include "bitmap.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <GL/gl.h>
/**
* Convenient way to group image parameters which tend to go together.
*/
typedef struct Image8u {
/** The actual pixel data of the image. */
Ipp8u* pixels;
/**
* Interval in bytes between consecutive rows. This may be different
* than the image width because ippiMalloc adds padding for better
* alignment.
*/
int step;
/**
* Region of interest (image size).
*/
IppiSize roi;
} Image8u;
/**
* Allocate a new image.
*/
static Image8u new_Image8u(int width, int height) {
int step;
Ipp8u* pixels = ippiMalloc_8u_C1(width, height, &step);
IppiSize roi = {width, height};
Image8u out = {pixels, step, roi};
return out;
}
/**
* Free memory associated with an image.
*/
static void free_Image8u(Image8u image) {
ippiFree(image.pixels);
}
/**
* Convenient way to group image parameters which tend to go together.
*/
typedef struct Image16s {
/** The actual pixel data of the image. */
Ipp16s* pixels;
/**
* Interval in pixels between consecutive rows. Multiply this by
* sizeof(Ipp16s) before using in ippi functions.
*/
int step;
IppiSize roi;
} Image16s;
/**
* Allocate a new image.
*/
static Image16s new_Image16s(int width, int height) {
int step;
Ipp16s* pixels = ippiMalloc_16s_C1(width, height, &step);
IppiSize roi = {width, height};
Image16s out = {pixels, step/sizeof(Ipp16s), roi};
return out;
}
/**
* Free memory associated with an image.
*/
static void free_Image16s(Image16s image) {
ippiFree(image.pixels);
}
/**
* Allocate a new Line object.
*/
static Line* new_Line() {
Line* out = malloc(sizeof(Line));
out->num_pixels = 0;
out->num_votes = 0;
}
void print_lines(Line* lines, int num_lines) {
int i;
printf("%d Lines:\n", num_lines);
for (i = 0; i < num_lines; ++i) {
Line* line = &lines[i];
printf("Line from (%f %f) to (%f %f)\n",
line->x1, line->y1,
line->x2, line->y2);
}
}
void free_lines(Line* lines) {
free(lines);
}
/**
* Part of a connected component of an image. Connected components are
* represented by tree loops (i.e. n-ary trees where the root points to itself)
* where each node is a Region.
*/
typedef struct Region {
/** Pointer to the parent region. */
struct Region* parent;
/**
* Rank is a measure of the complexity of a connected component.
* It is used when merging two connected components to decide which
* will be the root node.
*/
int rank;
/** Line associated with the region. */
Line* line;
/**
* Regions are allocated in a linked list for easy iteration, this
* points to the next region in the list.
*/
struct Region* next;
} Region;
/**
* Allocate a new region / connected component.
*/
static Region* new_Region() {
Region* out = malloc(sizeof(Region));
out->parent = out;
out->rank = 0;
out->line = NULL;
out->next = NULL;
return out;
}
/**
* Find the root Region of a connected component. As a side-effect this updates
* the tree so that all the Regions on the path to the root point directly to
* the root, speeding up future lookups.
*/
static Region* region_root(Region* c) {
if (c->parent == c) return c;
return c->parent = region_root(c->parent);
}
/**
* Merge two connected components by linking one's root to the other. The
* smaller component is always linked to the larger one, to avoid
* badly unbalanced trees.
*/
static void union_regions(Region* a, Region* b) {
if (a == b) return;
Region* aroot = region_root(a);
Region* broot = region_root(b);
if (aroot->rank > broot->rank) {
broot->parent = aroot;
} else if (aroot->rank < broot->rank) {
aroot->parent = broot;
} else if (aroot != broot) {
broot->parent = aroot;
++aroot->rank;
}
}
/**
* Map of pixels to regions. This mainly exists for compatibility with the AT
* macro.
*/
typedef struct ImageR {
/** Array of pointers to regions. */
Region** pixels;
/** Interval (in units of sizeof(void*)) between consecutive rows. */
int step;
/** Region of interest; size of the image. */
IppiSize roi;
} ImageR;
/**
* Allocate a new region map.
*/
static ImageR new_ImageR(int width, int height) {
IppiSize roi = {width, height};
ImageR out = {calloc(width*height, sizeof(Region*)), width, roi};
return out;
}
/**
* Free data associated with a region map.
*/
static void free_ImageR(ImageR image) {
free(image.pixels);
}
/**
* Return the element at coordinates (x,y) as measured from the top left in an
* Image8u-like struct.
*/
#define AT(image, x, y) (image).pixels[(y)*(image).step + (x)]
/**
* Print the contents of a region map for debugging purposes.
*/
static void print_ImageR(ImageR image) {
int x, y;
for (y = 0; y < image.roi.height; ++y) {
for (x = 0; x < image.roi.width; ++x) {
printf("%p ", AT(image, x, y));
}
printf("\n");
}
}
/**
* Print the contents of a grayscale image for debugging purposes.
*/
static void print_Image8u(Image8u image) {
int x, y;
for (y = 0; y < image.roi.height; ++y) {
for (x = 0; x < image.roi.width; ++x) {
printf("%d", AT(image, x, y)/26);
}
printf("\n");
}
}
/**
* Calculate 4-connected components of the image based on gradient orientation.
* See http://people.csail.mit.edu/rahimi/connected/ for more details.
*/
static ImageR line_support(Image8u gradient_m, Image8u gradient_q, Region** rtail) {
int x, y;
ImageR out = new_ImageR(gradient_m.roi.width, gradient_m.roi.height);
Region* rhead = *rtail;
// first pixel
if (AT(gradient_m, 0, 0)) {
AT(out, 0, 0) = (*rtail) = (*rtail)->next = new_Region();
}
// first row
for (x = 1; x < gradient_m.roi.width; ++x) {
if (AT(gradient_m, x, 0)) {
if (AT(out, x-1, 0) && AT(gradient_q,x,0) == AT(gradient_q,x-1,0)) {
// matches left pixel
AT(out, x, 0) = AT(out, x-1, 0);
} else {
AT(out, x, 0) = (*rtail) = (*rtail)->next = new_Region();
}
}
}
// remaining rows
for (y = 1; y < gradient_m.roi.height; ++y) {
// first pixel of row
if (AT(gradient_m, 0, y)) {
if (AT(out, 0, y-1) && AT(gradient_q,0,y) == AT(gradient_q,0,y-1)) {
// matches top pixel
AT(out, 0, y) = AT(out, 0, y-1);
} else {
AT(out, 0, y) = (*rtail) = (*rtail)->next = new_Region();
}
}
// remaining pixels of row
for (x = 1; x < gradient_m.roi.width; ++x) {
if (AT(gradient_m, x, y)) {
if (AT(out, x-1, y) && AT(gradient_q,x,y) == AT(gradient_q,x-1,y)) {
// matches left pixel
AT(out, x, y) = AT(out, x-1, y);
if (AT(out, x, y-1) && AT(gradient_q,x,y) == AT(gradient_q,x,y-1)) {
// matches top and left pixels
union_regions(AT(out, x, y), AT(out, x, y-1));
AT(out, x, y) = AT(out, x, y-1);
}
} else if (AT(out, x, y-1) && AT(gradient_q,x,y-1) == AT(gradient_q,x,y)) {
// matches top pixel
AT(out, x, y) = AT(out, x, y-1);
} else {
// matches neither top nor left pixels
AT(out, x, y) = (*rtail) = (*rtail)->next = new_Region();
}
}
}
}
// Merge regions and compute dimensions of each region
for (y = 0; y < gradient_m.roi.height; ++y) {
for (x = 0; x < gradient_m.roi.width; ++x) {
if (!AT(out, x, y)) continue;
// Find the root for this connected component.
AT(out, x, y) = region_root(AT(out, x, y));
Region* region = AT(out, x, y);
if (!region->line) {
region->line = new_Line();
region->line->x1 = region->line->x2 = x+1;
region->line->y1 = region->line->y2 = y+1;
} else {
if (x+1> region->line->x2) region->line->x2 = x+1;
if (x+1< region->line->x1) region->line->x1 = x+1;
region->line->y2 = y+1;
}
++region->line->num_pixels;
}
}
// remove temporary regions from the linked list
Region* rtmp;
while (rhead->next != NULL) {
if (!rhead->next->line) {
rtmp = rhead->next->next;
free(rhead->next);
rhead->next = rtmp;
} else {
rhead = rhead->next;
}
}
*rtail = rhead;
return out;
}
static void estimate_line_properties(ImageR regions, Image8u gradient_m, Image16s gradient_x, Image16s gradient_y) {
int x, y;
// calculate correct line parameters for good lines
for (y = 0; y < regions.roi.height; ++y) {
for (x = 0; x < regions.roi.width; ++x) {
Region* region = AT(regions, x, y);
if (!region) continue;
Line* line = region->line;
if (!line) continue;
int new_num_pixels = line->num_pixels + 1;
double weight = hypot(line->tangent_x, line->tangent_y) * line->num_pixels;
double sx = AT(gradient_x, x, y);
double sy = AT(gradient_y, x, y);
line->tangent_x = line->tangent_x * line->num_pixels / (double)new_num_pixels + sx / new_num_pixels;
line->tangent_y = line->tangent_y * line->num_pixels / (double)new_num_pixels + sy / new_num_pixels;
double new_weight = hypot(line->tangent_x, line->tangent_y) * new_num_pixels;
double current_weight = hypot(sx, sy);
line->midpoint_x = (line->midpoint_x * weight + (x+1) * current_weight) / (weight + current_weight);
line->midpoint_y = (line->midpoint_y * weight + (y+1) * current_weight) / (weight + current_weight);
line->num_pixels = new_num_pixels;
}
}
}
static void draw_Image8u(Image8u image) {
glClear(GL_COLOR_BUFFER_BIT);
glColor3f(1,0,0);
glLineWidth(2);
glRasterPos2f(0,0);
glDrawPixels(image.roi.width, image.roi.height, GL_LUMINANCE, GL_UNSIGNED_BYTE, image.pixels);
}
static int min(int a, int b) {
if (a <= b) return a;
else return b;
}
static int max(int a, int b) {
if (a >= b) return a;
else return b;
}
void burns_line_extraction(Ipp8u* pixels, int step, int width, int height,
int num_buckets, int min_gradient, double min_length,
Line** lines, int* num_lines)
{
IppiSize roi = {width, height};
Image8u image = {pixels, step, roi};
int x, y;
// The edges of the image do not have a well-defined gradient (using
// the Sobel filter), so we will ignore them for the purposes of this
// algorithm.
roi.height -= 2;
roi.width -= 2;
Image16s gradient_x = new_Image16s(roi.width, roi.height);
ippiFilterSobelHoriz_8u16s_C1R(&AT(image, 1, 1), image.step, gradient_x.pixels, gradient_x.step*sizeof(Ipp16s), gradient_x.roi, ippMskSize3x3);
Image16s gradient_y = new_Image16s(roi.width, roi.height);
ippiFilterSobelVert_8u16s_C1R(&AT(image, 1, 1), image.step, gradient_y.pixels, gradient_y.step*sizeof(Ipp16s), gradient_y.roi, ippMskSize3x3);
// Gradient magnitude and direction (in integer multiples of pi/8).
// We could do this with ippiFilter but since we have to iterate
// through the image anyways, doing the filtering at the same time
// avoids unnecessary memory allocation.
Image8u gradient_m = new_Image8u(roi.width, roi.height);
// Gradient direction quantized into buckets starting at 0 degrees.
Image8u gradient_q1 = new_Image8u(roi.width, roi.height);
// Gradient direction quantized into buckets starting at 180/num_buckets
// degrees (i.e. offset by half a bucket from gradient_q1).
Image8u gradient_q2 = new_Image8u(roi.width, roi.height);
for (y = 0; y < gradient_m.roi.height; ++y) {
for (x = 0; x < gradient_m.roi.width; ++x) {
int hgrad = AT(gradient_x, x, y);
int vgrad = AT(gradient_y, x, y);
int gmag = abs(hgrad)/4+ abs(vgrad)/4;
if (gmag > min_gradient) {
double theta = atan2((double)vgrad, (double)hgrad);
int bucket = (int)floor((theta+M_PI) * (double)num_buckets/M_PI);
AT(gradient_m, x, y) = gmag;
AT(gradient_q1, x, y) = bucket/2;
AT(gradient_q2, x, y) = ((bucket+1)%(2*num_buckets))/2;
} else {
AT(gradient_m, x, y) = 0;
AT(gradient_q1, x, y) = 0;
AT(gradient_q2, x, y) = 0;
}
}
}
// Find the line-support regions of the image; regions will be
// collected in a linked list so we can iterate through them later.
Region* rhead = new_Region();
Region* rtail = rhead;
ImageR regions1 = line_support(gradient_m, gradient_q1, &rtail);
ImageR regions2 = line_support(gradient_m, gradient_q2, &rtail);
// free temporary list head
rtail = rhead;
rhead = rhead->next;
free(rtail);
// gradient orientations no longer needed
free_Image8u(gradient_q1);
free_Image8u(gradient_q2);
// approximate length of each region's line
for (rtail = rhead; rtail != NULL; rtail = rtail->next) {
Line* line = rtail->line;
line->length = abs(line->x1 - line->x2)/2+ abs(line->y1 - line->y2)/2;
}
// vote on regions
for (y = 0; y < regions1.roi.height; ++y) {
for (x = 0; x < regions1.roi.width; ++x) {
if (!AT(regions1, x, y)) continue;
Region* c1 = AT(regions1, x, y);
Region* c2 = AT(regions2, x, y);
if (c1->line->length >= c2->line->length) {
++c1->line->num_votes;
} else {
++c2->line->num_votes;
}
}
}
// discard bad (undervoted) lines
for (rtail = rhead; rtail != NULL; rtail = rtail->next) {
Line* line = rtail->line;
if (line) {
if (line->length < min_length || line->num_votes < line->num_pixels/2) {
free(line);
rtail->line = NULL;
} else {
// reset pixel counter for averages
line->num_pixels = 0;
}
}
}
estimate_line_properties(regions1, gradient_m, gradient_x, gradient_y);
estimate_line_properties(regions2, gradient_m, gradient_x, gradient_y);
// calculate more accurate endpoints for lines
for (rtail = rhead; rtail != NULL; rtail = rtail->next) {
Line* line = rtail->line;
if (line) {
double xmin = line->x1; double xmax = line->x2;
double ymin = line->y1; double ymax = line->y2;
line->x1 = line->midpoint_x - (line->midpoint_y - ymin)*line->tangent_x/line->tangent_y;
if (line->x1 < xmin) {
line->x1 = xmin;
line->y1 = line->midpoint_y - (line->midpoint_x - xmin)*line->tangent_y/line->tangent_x;
} else if (line->x1 > xmax) {
line->x1 = xmax;
line->y1 = line->midpoint_y + (xmax - line->midpoint_x)*line->tangent_y/line->tangent_x;
}
line->x2 = line->midpoint_x + (ymax - line->midpoint_y)*line->tangent_x/line->tangent_y;
if (line->x2 < xmin) {
line->x2 = xmin;
line->y2 = line->midpoint_y - (line->midpoint_x - xmin)*line->tangent_y/line->tangent_x;
} else if (line->x2 > xmax) {
line->x2 = xmax;
line->y2 = line->midpoint_y + (xmax - line->midpoint_x)*line->tangent_y/line->tangent_x;
}
// re-calculate length
line->length = hypot(line->x1 - line->x2, line->y1 - line->y2);
// discard too-short lines
if (line->length < min_length) {
free(line);
rtail->line = NULL;
}
}
}
// count lines
*num_lines = 0;
for (rtail = rhead; rtail != NULL; rtail = rtail->next) {
Line* line = rtail->line;
if (line) ++(*num_lines);
}
// gradient information no longer needed
free_Image8u(gradient_m);
free_Image16s(gradient_y);
free_Image16s(gradient_x);
// region maps no longer needed
free_ImageR(regions1);
free_ImageR(regions2);
// Lines will be collected in an array
*lines = malloc(*num_lines * sizeof(Line));
int i = 0;
// build output list and free regions
while (rhead != NULL) {
Line* line = rhead->line;
if (line) {
// found the root of a connected component with a good line
(*lines)[i++] = *line;
free(line);
}
// free the region and move to the next one
rtail = rhead;
rhead = rhead->next;
free(rtail);
}
}
void burns_line_extraction_demo(Line** lines, int* num_lines) {
IppiSize roi;
int rgb_step, gray_step;
Ipp8u* rgb_pixels = read_bmp("demo.bmp", &roi.width, &roi.height, &rgb_step);
Ipp8u* gray_pixels = ippiMalloc_8u_C1(roi.width, roi.height, &gray_step);
ippiRGBToGray_8u_C3C1R(rgb_pixels, rgb_step, gray_pixels, gray_step, roi);
burns_line_extraction(gray_pixels, gray_step, roi.width, roi.height, 8, 50, 10, lines, num_lines);
}
最后是调用和利用OpenGL进行显示的主函数main.c
#include "burns.h"
#include "bitmap.h"
#include <stdio.h>
#include <math.h>
#include <ippcc.h>
#include <GL/gl.h>
#include <GL/glut.h>
Ipp8u* rgb_pixels;
int rgb_step;
IppiSize rgb_roi;
Ipp8u* gray_pixels;
int gray_step;
IppiSize gray_roi;
void display() {
glClear(GL_COLOR_BUFFER_BIT);
glLineWidth(2);
glRasterPos2f(0,0);
glDrawPixels(gray_roi.width, gray_roi.height, GL_LUMINANCE, GL_UNSIGNED_BYTE, gray_pixels);
Line* lines;
int num_lines;
int i;
burns_line_extraction(gray_pixels, gray_step, gray_roi.width, gray_roi.height,
8, 10, 5,
&lines, &num_lines);
glColor3f(1,0,0);
glBegin(GL_LINES);
for (i = 0; i < num_lines; ++i) {
Line* line = &lines[i];
glVertex3d(line->x1, line->y1, 0);
glVertex3d(line->x2, line->y2, 0);
}
glEnd();
glPointSize(3);
glColor3f(0,0,1);
glBegin(GL_POINTS);
for (i = 0; i < num_lines; ++i) {
Line* line = &lines[i];
glVertex3d(line->midpoint_x, line->midpoint_y, 0);
}
glEnd();
glFlush();
free_lines(lines);
}
int main(int argc, char **argv) {
IppiSize rgb_roi;
rgb_pixels = read_bmp("demo.bmp", &rgb_roi.width, &rgb_roi.height, &rgb_step);
gray_pixels = ippiMalloc_8u_C1(rgb_roi.width, rgb_roi.height, &gray_step);
gray_roi = rgb_roi;
ippiRGBToGray_8u_C3C1R(rgb_pixels, rgb_step, gray_pixels, gray_step, rgb_roi);
glutInit(&argc, argv);
glutInitWindowSize(gray_roi.width, gray_roi.height);
glutInitWindowPosition(0, 0);
glutInitDisplayMode(GLUT_RGB);
glutCreateWindow("window");
glClearColor(0.0, 0.0, 0.0, 0.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(0, gray_roi.width, 0, gray_roi.height, -1, 1);
glutDisplayFunc(display);
glutMainLoop();
return 0;
}
编译过程中可能遇到的问题
理论上这几个c文件利用gcc是可以直接编译通过的，但无奈本人的gcc功力不够，最后只能把它移植到VC6中进行编译和调试。下面列 一下用VC6编译过程中可能遇到的问题和解决方案：
1、在自定义的min和max函数报错：错误：expected identifier or ‘(’ before ‘int’，找了半天感觉应该是max或者min这个函数名称已经被用过了？改成Mymax或者Mymin以后就行了
2、结构体初始化报错，比如Image8u out = {pixels, step, roi};改成各个成员分别初始化就行了：
Image8u out;
out.pixels = pixels;
out.step = step;
out.roi = roi;
3、OpenGL环境的配置可以参考我之前写过的博客Windows下配置OpenGL的开发环境，以VC6为例
需要的资源下载
最后列一下编译需要的资源，其实主要是ipp库和OpenGL
OpenGL上面说过了，参考我之前发过的博客，里面有资源下载地址，找个简单例子测试下没问题就行了。
ipp库就比较恶心了，我是在intel的网站下载了个试用版http://software.intel.com/en-us/intel-ipp
反正有头文件和lib就行了。
最后贴几张效果图吧，费这么大劲总得搞点东西出来啊。
需要注意的是不知为啥最后只能显示宽度为64倍数的bmp图片（其他宽度的图片提取结果应该是没问题，只是显示一团糟，怀疑是bitmap的数据结构导致的，高手可以改进下，争取去掉它的ipp数据结构，换成opencv或者gdal的数据结构，这样就啥图片都可以处理了）
首先是经典的lena中的直线
再来看一张高分辨率遥感数据的直线提取结果
