How to determine if a list of polygon points are in clockwise order?
Some of the suggested methods will fail in the case of a non-convex polygon, such as a crescent. Here's a simple one that will work with non-convex polygons (it'll even work with a self-intersecting polygon like a figure-eight, telling you whether it's mostly clockwise).
Sum over the edges, (x2 − x1)(y2 + y1). If the result is positive the curve is clockwise, if it's negative the curve is counter-clockwise. (The result is twice the enclosed area, with a +/- convention.)
The intersection of the planes can be non-robust when the planes are near parallel. Here are the tricks we use in our implementation:
We first test if an intersection exists for every pair of planes. Given N planes, a naive pairwise intersection algorithm has a complexity of O(N 2 ). We then collect plane triplets such that every pair in the plane triplet intersect. This is achieved by testing each plane against the known intersecting pairs, which also has a complexity of O(N2 ).
The 3D vertices of the final faces are obtained by computing the intersections of the plane triplets. To cope with limited floating point precision, each vertex is identified by the indices (in an increasing order) of the three planes from which it is computed. By doing so, two vertices with almost identical positions can be distinguished. This turned out to be quite robust in handling very close and near parallel planes.
For more details, please refer to our implementation (i.e., hypothesis_generator.h and hypothesis_generator.cpp).
For your own data, you can use my Mapple to extract planes. Here is the link to Mapple:https://3d.bk.tudelft.nl/liangliang/software.html After you load the point cloud to Mapple, go to 'Partition' menu and click 'Extract Primitives'. To visualize the planes, change the renderer from 'Plain' to 'Group' in the Rendering panel (at the left side of Mapple). You can save the planes as bvg (Binary Vertex Group) format. The ASCII format vg also works but slow.
Below you will find the description of the file format and a simple example file.
File format:
PolyFit assumes that planar segments have been extracted properly and are stored in vg (vertex group) format I have developed for my research projects. The general vg format allows you to save a point cloud followed by its segmentation information:
For each point, its coordinates, normal, and color.
For each segment (a group of points representing a primitive or an object), its label, model parameters, color, and the indices of the points that belong to this segment.
PolyFit handles planes only, so the model parameters are simply the plane parameters (e.g., a, b, c, and d in the plane equation ax + by + cz + d = 0).
Below is a details description of the ASCII vg format. The source code of PolyFit also contains an implementation for binary format. Please refer to 'point_set_serializer_vg.cpp' of the source code for more information.
// the coordinates of the points
num_points: N // N is an integer denoting the number of points
x1 y1 z1 // 3 floating point numbers
...
xN yN zN
// the colors of the points
num_colors: N // N is an integer denoting the number of colors (can be different from the number of points)
r1 g1 b1 // 3 floating point numbers
...
rN gN bN
// the normals of the points
num_normals: N // N is an integer denoting the number of normals (can be different from the number of points or colors)
nx1 ny1 nz1 // 3 floating point numbers
...
nxN nyN nzN
// now we store the segmentation information
num_groups: M // M is an integer denoting the number of segments/primitives/objects in this point cloud
// now the information for the 1st segment/primitive/object
group_type: 0 // must be 0 standing for a plane (1 to 5 are other types of primitives)
num_group_parameters: 4 // must be 4 (planes are represented using 4 parameters)
group_parameters: a b c d // 4 floating point numbers (e.g., a, b, c, and d for a plane)
group_label: label // the label (a string) of the 1st vertex group
group_color: r g b // 3 floating point numbers denoting the color of the 1st vertex group
group_num_points: N // N is an integer denoting the number of points in the 1st vertex group
id1 ... idN // N integer numbers denoting the indices of the points in the 1st vertex group
num_children: 0 // a segment/primitive/object may contain subgroups, but for PolyFit this must be 0
...
group_type: 0 // here starts the last segment in the point cloud (similar to the 1st one)
num_group_parameters: 4
group_parameters: a b c d
group_label: label
group_color: r g b
group_num_points: N
id1 ... idN
num_children: 0
A tiny example file (comments are NOT part of the file):
num_points: 210882 // there are 210882 points in this point cloud
-0.06324 0.03597 0.04208
...
-0.06449 0.03651 0.04043
num_colors: 210882
1 0.04 0
...
0 0.91 0
num_normals: 210882
-0.65573 -0.50319 0.56283
...
-0.56256 -0.51472 0.64698
num_groups: 7 // there are 7 segments in this point cloud
group_type: 0 // the first segment is a plane
num_group_parameters: 4 // there are 4 parameters for the plane
group_parameters: -0.00494 -0.11430 0.99343 -0.03321 //a,b,c,and d in the plane equation ax+by+cz+d=0
group_label: unknown // the label for the plane is ¡°unknown¡± (not given)
group_color: 0.1 0.6 0.2
group_num_point: 14433 // the plane consists of 14433 points
30493 ... 8798 // here are the indices of the 14433 points
num_children: 0
... // here are 5 segments/primitives/objects
group_type: 0 // here comes the last segment. It is also a plane.
num_group_parameters: 4
group_parameters: 0.67223 0.25087 0.69654 -0.02011
group_label: unknown
group_color: 0.7 0.2 0.4
group_num_point: 9155
38812 ... 140417
num_children: 0
Parameters.
The parameters for most examples are as follows: fitting = 0.46, coverage = 0.27, and complexity = 0.27 (Note that weights in a wide range can produce the same results). Slightly different weights (fitting = 0.3, coverage = 0.4, and complexity = 0.3) are used for the sofa example in Figure 4(i), where the background (ground plane) has a much higher density than the object (sofa), thus the smaller data fitting weight. In case non-default parameters are used, these parameters are provided in the files names.
- 모듈들이 특정한 실행 파일 안에 함께 포함되어 있다면 이것은 명확히 하나의 프로그램으로 결합되어 있는 것
- 파이프와 소켓, 명령행 인자 등은 두개의 독립된 프로그램간의 통신을 위해서 사용되는 매커니즘입니다. 따라서 모듈들이 이러한 형식을 사용한다면 모듈들은 독립된 프로그램으로 볼 수 있습니다.
GPL 라이센스가 걸린 소스의 결과물과 링킹을 했다고 소스를 공개해야 하는것은 아닙니다. (위에서는 무조건 공개해야 하는것 처럼 적혀 있지요.)
Linux의 core들도 GPL로 되어 있다고 본것 같은데, 그렇다면 Linux core(API)를 사용하는것이면 모두 공개 되어야 하는게 아닌가? 라고 생각할수 있는데 OS에 밀접하게 연관되어 있는 (main component?)에 대해서는 문제가 되지 않는다고 합니다.
OpenCV의 cvCanny() 함수는 이러한 과정 없이 원본 이미지에서 바로 canny edge를 찾는다. OpenCV를 사용할 때 노이즈에 의한 잡음이 많은 경우 Gaussian blur를 적용하고 cvCanny()를 사용해야겠다. 아무래도 OpenCV에서 Gaussian blur를 적용하지 않는것은 연산량이 많기 때문인것 같은데 Separable Kernel을 사용하면 속도 향상에 도움이 된다.Separable Kernel을 이용한 Gaussian blur에 대해서는 아래 주소에 포스팅에 정리해 두었다.
엣지의 magnitude 가 maximum인 경우만 엣지로 사용한다. 이 때 모든 3x3 neighbor에 대해서 maximum인지 판별하는 것이 아닌 edge의 orientaton을 고려해 탐색 방향을 결정한다. 엣지 방향은 아래 그림과 같이 4개로 나누며
maximum인지 판단하는 방향은 방향이 0이면 좌우, 1이면 현재 위치에서 offset이 (-1,-1), (+1,+1)인 점과 비교한다. 방향이 2라면 위, 아래와 비교, 3이면 offset이 (-1,+1), (+1,-1)인 점과 비교한다.
Canny Edge의 특징중 하나는 두개의 threshold를 사용하는 것인데, 이 과정에서 높은 threshold보다 엣지의 magnitude가 크면 엣지 픽셀로 지정하고 낮은 threshold 보다 높은 magnitude를 가지면 나중에 edge들을 연결(hysteresis analysis)할 때 엣지 픽셀과 인접해 있으면 엣지로 지정한다.
다음 그림은 non-maxima suppression을 수행한 결과로 maxima는 흰색으로, non-maxima는 회색으로, 엣지가 아닌 부분은 검은색으로 표시했다. 이 그림은 단지 중간 과정을 보여주기 위한 것으로 non-maxima 픽셀도 역시 엣지픽셀로 사용하지 않는다. 또한 threshold를 적용하지 않고 단지 엣지픽셀들의 orientation에 따라 magnitude가 가장 큰 픽셀들만 보여주는 것이다.
4. 엣지 연결하기
위에서도 말했듯이 Canny edge detector는 두개의 threshold를 사용한다. 그 이유는 큰 threshold 보다 큰 magnitude를 가지는 엣지들만 엣지로 할 경우 영상의 잡음으로 인해 엣지 픽셀이 엣지픽셀로 인식되지 않을 경우가 있어 낮은 threshold 보다 높은 magnitude를 가지고 엣지 픽셀과 인접한 픽셀들을 잇기 위함이다. 편의상 큰 threshold를 T1, 작은 threshold를 T2라고 부르겠다.
다음 그림에서 왼쪽은 non-maximum suppression 후 T1 보다 높은 edge픽셀만을 보여주는 것이고 오른쪽은 T2보다 magnitude가 크고 엣지와 인접한 픽셀들을 연결한 결과이다. 이 때 T1=100, T2=50을 사용했다.
이 과정에서는 모든 엣지에 대해서 주변에 인접하고 magnitude가 T2보다 높은 픽셀들을 엣지로 지정하고 또 새롭게 엣지로 지정된 픽셀 주변에서 magnitude가 T2보다 높은 인접한 픽셀을 엣지로 지정하는 식으로 반복적으로 수행한다. OpenCV의 구현에서는 이 부분에 stack을 사용했다. 이 stack의 pop은 말 그대로 pop인데 push는 push할 픽셀을 엣지로 지정하고 난 후 stack에 저장한다. 내가 구현한 코드에서도 이 부분을 따라했다.
상세한 설명을 하자면 위 3번 단계에서 T1보다 큰 픽셀은 픽셀로 지정하고 이 픽셀의 주소를 stack에 push한다. 그리고 T1보다 작고 T2보다 큰 픽셀은 그냥 엣지가 될 수 있는 것임을 지정한다. 이 과정이 모두 끝나면 stack에서 엣지 픽셀을 꺼내고 해당 엣지 픽셀의 3x3 neighbor에 T2보다 큰 픽셀이 있을 경우 이 픽셀을 엣지로 지정하고 stack에 push한다. 이 과정을 stack이 비게 될 때 까지 반복하면 엣지가 모두 연결이 된다.
그리고 gaussian blur 를 적용한 것과 그렇지 않은 것의 차이를 비교하기 위해 gaussian blur를 적용하지 않고 canny edge를 찾은 결과를 보면 다음과 같다.
다음은 내가 C로 구현한 canny edge detector source code이다. 함수 인자 중 이미지와 결과 이미지는 BYTE* 로 주고받는데 OpenCV를 사용하는 경우라면 IplImage의 imageData 를 사용하면 된다.
그런데 어디에서 속도차이가 나는지는 모르겠지만 OpenCV의 코드가 내 코드보다 0.5ms 정도 더 빠르다. 내 코드는 320x240크기의 이미지에서 대략 2.5ms 정도 걸리지만 OpenCV는 2ms 정도 걸리는데 아무리 찾아봐도 무엇때문에 이런 속도 차이가 나는지 모르겠다. 분명 Intel 의 CPU 구조를 잘 아는 사람들이 만든거니 더 최적화된 코드를 짰을텐데 말이다. 그래서 그런지 Intel의 C 컴파일러와 MS의 C컴파일러로 같은 코드를 컴파일해도 Intel의 C컴파일러가 훨씬 더 빠른 최적화된 코드를 만든다고 하는데 Intel C 컴파일러는 무료로 제공되지 않아 사용해보질 못해 아쉽다. MS의 dreamspark같이 학생들에게 무료로 컴파일러를 제공하면 좋을텐데...
void CannyEdge(BYTE *pImage, int width, int height, int th_high, int th_low, BYTE *pEdge)
{
int i, j;
int dx, dy, mag, slope, direction;
int index, index2;
const int fbit = 10;
const int tan225 = 424; // tan25.5 << fbit, 0.4142
const int tan675 = 2472; // tan67.5 << fbit, 2.4142
const int CERTAIN_EDGE = 255;
const int PROBABLE_EDGE = 100;
bool bMaxima;
int *mag_tbl = new int[width*height];
int *dx_tbl = new int[width*height];
int *dy_tbl = new int[width*height];
BYTE **stack_top, **stack_bottom;
stack_top = new BYTE*[width*height];
stack_bottom = stack_top;
#define CANNY_PUSH(p) *(p) = CERTAIN_EDGE, *(stack_top++) = (p)
#define CANNY_POP() *(--stack_top)
for(i=0; i<width*height; i++) {
mag_tbl[i] = 0;
pEdge[i] = 0;
}
// Sobel Edge Detection
for(i=1; i<height-1; i++) {
index = i*width;
for(j=1; j<width-1; j++) {
index2 = index+j;
// -1 0 1
// -2 0 2
// -1 0 1
dx = pImage[index2-width+1] + (pImage[index2+1]<<1) + pImage[index2+width+1]
-pImage[index2-width-1] - (pImage[index2-1]<<1) - pImage[index2+width-1];
// -1 -2 -1
// 0 0 0
// 1 2 1
dy = -pImage[index2-width-1] - (pImage[index2-width]<<1) - pImage[index2-width+1]
+pImage[index2+width-1] + (pImage[index2+width]<<1) + pImage[index2+width+1];
mag = abs(dx)+abs(dy); // magnitude
//mag = sqrtf(dx*dx + dy*dy);
dx_tbl[index2] = dx;
dy_tbl[index2] = dy;
mag_tbl[index2] = mag; // store the magnitude in the table
} // for(j)
} // for(i)
for(i=1; i<height-1; i++) {
index = i*width;
for(j=1; j<width-1; j++) {
index2 = index+j;
mag = mag_tbl[index2]; // retrieve the magnitude from the table
// if the magnitude is greater than the lower threshold
if(mag > th_low) {
// determine the orientation of the edge
dx = dx_tbl[index2];
dy = dy_tbl[index2];
if(dx != 0) {
slope = (dy<<fbit)/dx;
if(slope > 0) {
if(slope < tan225)
direction = 0;
else if(slope < tan675)
direction = 1;
else
direction = 2;
}
else {
if(-slope > tan675)
direction = 2;
else if(-slope > tan225)
direction = 3;
else
direction = 0;
}
}
else
direction = 2;
bMaxima = true;
// perform non-maxima suppression
if(direction == 0) {
if(mag < mag_tbl[index2-1] || mag < mag_tbl[index2+1])
bMaxima = false;
}
else if(direction == 1) {
if(mag < mag_tbl[index2+width+1] || mag < mag_tbl[index2-width-1])
bMaxima = false;
}
else if(direction == 2){
if(mag < mag_tbl[index2+width] || mag < mag_tbl[index2-width])
bMaxima = false;
}
else { // if(direction == 3)
if(mag < mag_tbl[index2+width-1] || mag < mag_tbl[index2-width+1])
bMaxima = false;
}
if(bMaxima) {
if(mag > th_high) {
pEdge[index2] = CERTAIN_EDGE; // the pixel does belong to an edge
*(stack_top++) = (BYTE*)(pEdge+index2);
}
else
pEdge[index2] = PROBABLE_EDGE; // the pixel might belong to an edge
}
}
} // for(j)
} // for(i)
while(stack_top != stack_bottom) {
BYTE* p = CANNY_POP();
if(*(p+1) == PROBABLE_EDGE)
CANNY_PUSH(p+1);
if(*(p-1) == PROBABLE_EDGE)
CANNY_PUSH(p-1);
if(*(p+width) == PROBABLE_EDGE)
CANNY_PUSH(p+width);
if(*(p-width) == PROBABLE_EDGE)
CANNY_PUSH(p-width);
if(*(p-width-1) == PROBABLE_EDGE)
CANNY_PUSH(p-width-1);
if(*(p-width+1) == PROBABLE_EDGE)
CANNY_PUSH(p-width+1);
if(*(p+width-1) == PROBABLE_EDGE)
CANNY_PUSH(p+width-1);
if(*(p+width+1) == PROBABLE_EDGE)
CANNY_PUSH(p+width+1);
}
for(i=0; i<width*height; i++)
if(pEdge[i]!=CERTAIN_EDGE)
pEdge[i] = 0;
delete [] mag_tbl;
delete [] dx_tbl;
delete [] dy_tbl;
delete [] stack_bottom;
}
PolygonFill() 에서 폴리곤의 최대 최소 범위의 픽셀들을 하나씩 체크하면서 내부, 외부를 판단함.
내부 외부 판단은 floodFill4() 를 이용함.
// ***************************************
// CS459-Fundamentals of Computer Graphics
// Project 1: Interactive Polygon-Filling
// Student : Ahmad Pahlavan Tafti
// ***************************************
#include <stdlib.h>
#include <GL/glut.h>
#include <math.h>
#include <windows.h>
#include <stdio.h>
class scrPt { // Declaring the points for polygon
public:
GLint x, y;
};
#define width 450 // Window Size
#define height 450 // Window Size
#define PI 3.14159 // PI Number
#define TWOPI PI*2 // 2PI
//n=number of poly
#define n 4
scrPt poly[n];
float currentBuffer[width][height][3];
void init (void)
{
glClearColor (1.0, 1.0, 1.0, 0.0); // Set display-window color to white.
glMatrixMode (GL_PROJECTION); // Set projection parameters.
gluOrtho2D (0.0, 200.0, 0.0, 150.0);
}
//calculate angle of two point to determine is inside polygon or not:
/* Philippe Reverdy presented a solution to compute the sum
of the angles made between the test point and each pair of points making
up the polygon. If this sum is 2pi then the point is an interior point,
if 0 then the point is an exterior point. This also works for polygons
with holes given the polygon is defined with a path made up of coincident
edges into and out of the hole as is common practice in many CAD packages.
Reference: www.eecs.umich.edu/eecs/courses/eecs380/HANDOUTS/PROJ2/InsidePoly.html */
double Angle2D(double x1, double y1, double x2, double y2)
{
double dtheta,theta1,theta2;
theta1 = atan2(y1,x1);
theta2 = atan2(y2,x2);
dtheta = theta2 - theta1;
while (dtheta > PI)
dtheta -= TWOPI;
while (dtheta < -PI)
dtheta += TWOPI;
return(dtheta);
}
//mask1=R(red) of mask table
float mask1[3][3]={
{1,0,1},
{1,0,1},
{1,1,1}};
//mask1=G(green) of mask table
float mask2[3][3]={
{0,1,1},
{1,0,1},
{1,1,1}};
//mask1=B(blue) of mask table
float mask3[3][3]={
{0,0,1},
{0,1,1},
{1,1,1}};
//get the color of pixel according to the mask table
void getColor(int x, int y, float drawColor[3])
{
drawColor[0]=mask1[x/30-(3*(x/90))][y/30-(3*(y/90))];
drawColor[1]=mask2[x/30-(3*(x/90))][y/30-(3*(y/90))];
drawColor[2]=mask3[x/30-(3*(x/90))][y/30-(3*(y/90))];
}
//check a point by Angle2d function whether is in polygon or not. *polygon is Vertices of polygon
int floodFill4 (scrPt *polygon,scrPt p)
{
int i;
double angle=0;
scrPt p1,p2;
for (i=0;i<n;i++) {
p1.x = polygon[i].x - p.x;
p1.y = polygon[i].y - p.y;
p2.x = polygon[(i+1)%n].x - p.x;
p2.y = polygon[(i+1)%n].y - p.y;
angle += Angle2D(p1.x,p1.y,p2.x,p2.y);
}
if (abs(angle) < PI)
return(FALSE);
else
return(TRUE);
}
int minx,miny,maxx,maxy;
//Draw function:
void PolygonFill (void)
{
int i,j;
glClear (GL_COLOR_BUFFER_BIT); // Clear display window.
glColor3f (1.0, 0.0, 0.0);
scrPt p1;
//First of all:determine the rectangle that polygon inserted in it:
//purpose: avoiding calculate additional points that we sure dont lie in polygon area.
minx=width+1;
for(i=0;i<n;i++)
if(poly[i].x<minx)minx=poly[i].x;
miny=height+1;
for(i=0;i<n;i++)
if(poly[i].y<miny)miny=poly[i].y;
maxx=-1;
for(i=0;i<n;i++)
if(poly[i].x>maxx)maxx=poly[i].x;
maxy=-1;
for(i=0;i<n;i++)
if(poly[i].y>maxy)maxy=poly[i].y;
//Second: check pixels of polygon rectangle
for(i=0;i<width;i++)
for(j=0;j<height;j++)
{
p1.x=i;
p1.y=j;
if((p1.x<=maxx&&p1.x>=minx)&&(p1.y>=miny&&p1.y<=maxy))
{
if(floodFill4(poly,p1))
{
//if true should determine it's color according to the mask table by getcolor function
getColor(p1.x,p1.y,currentBuffer[i][j]);
}
else
{
//else it is black
currentBuffer[i][j][0]=0;
currentBuffer[i][j][1]=0;
currentBuffer[i][j][2]=0;
}
}
else
{
//else not in polygon rect is black:
currentBuffer[i][j][0]=0;
currentBuffer[i][j][1]=0;
currentBuffer[i][j][2]=0;
}
}
//draw buffer on display
glDrawPixels(width, height, GL_RGB, GL_FLOAT, currentBuffer);
glFlush ( ); // Process all OpenGL routines as quickly as possible.
}
void spec ( int key, int x, int y )
{
int i;
switch ( key )
{
case GLUT_KEY_UP :
printf ("UP arrow pressed\n");
for(i=0;i<n;i++)
poly[i].x+=10;
PolygonFill();
break;
case GLUT_KEY_DOWN :
printf ("DOWN arrow pressed\n");
for(i=0;i<n;i++)
poly[i].x-=10;
PolygonFill();
break;
case GLUT_KEY_RIGHT :
printf ("RIGHT arrow pressed\n");
for(i=0;i<n;i++)
poly[i].y+=10;
PolygonFill();
break;
case GLUT_KEY_LEFT :
printf ("LEFT arrow pressed\n");
for(i=0;i<n;i++)
poly[i].y-=10;
PolygonFill();
break;
}
}
void winReshapeFcn (int newWidth, int newHeight)
{
glMatrixMode (GL_PROJECTION);
glLoadIdentity ( );
gluOrtho2D (0.0, (GLdouble) newWidth, 0.0, (GLdouble) newHeight);
glClear (GL_COLOR_BUFFER_BIT);
}
void main (int argc, char** argv)
{
glutInit (&argc, argv); // Initialize GLUT.
glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB); // Set display mode.
glutInitWindowPosition (100, 100); // Set top-left display-window position.
glutInitWindowSize (width, height); // Set display-window width and height.
glutCreateWindow ("Project 1: Interactive Polygon-Filling"); // Create display window.
init ( ); // Execute initialization procedure.
poly[0].x = 200;
poly[0].y = 400;
poly[1].x = 100;
poly[1].y = 300;
poly[2].x = 100;
poly[2].y = 200;
poly[3].x = 200;
poly[3].y = 100;
poly[4].x = 300;
poly[4].y = 200;
poly[5].x = 300;
poly[5].y = 300;
glutDisplayFunc (PolygonFill); // Send graphics to display window.
glutReshapeFunc (winReshapeFcn);
glutSpecialFunc( spec ) ;
glutMainLoop ( ); // Display everything and wait.
}
가끔 함수 내에서만 사용하는 구조체를 선언하여 정렬하고 싶은데, 이를 위해선 구조체의 비교함수를 함수 내에서 선언해야 한다.
어렵게 찾은 방법이라 기록을 남겨둔다.
int main() {
// This declares a lambda, which can be called just like a function
auto print_message = [](std::string message)
{
std::cout << message << "\n";
};
// Prints "Hello!" 10 times
for(int i = 0; i < 10; i++) {
print_message("Hello!");
}
}
