ขั้นตอนวิธีแบบกรวย

ขั้นตอนวิธีแบบกรวย (Cone Algorithm) เป็นขั้นตอนวิธีที่รวดเร็ว ซึ่งสามารถระบุได้อย่างแม่นยำว่าอนุภาคใดบ้างที่อยู่บนพื้นผิวที่มีลักษณะเป็นก้อนเล็กๆ ใน 3 มิติ จากการคำนวณโดยทั่วไปทางเรขาคณิต แสดงให้เห็นว่าขั้นตอนวิธีแบบนี้มีประโยชน์อย่างยิ่งกับวิทยาศาสตร์ที่ที่มีคำนวณหาพื้นที่ผิว และวิทยาศาสตร์ที่ศึกษาเกี่ยวกับเรื่องของนาโนเทคโนโลยี
== แนวคิดพื้นฐาน ==

สำหรับกลุ่มก้อนสามมิติที่ให้มาในรูปด้านขวามือ ซึ่งประกอบด้วยชุดของอนุภาค เราจะบอกได้อย่างไรว่าอนุภาคใดบ้างอยู่บนพื้นผิว? เหมือนจะง่ายถ้ามองจากตาของมนุษย์ เกณฑ์การมองเห็นอาจถูกอธิบายได้ เช่นถ้าเรากำหนดว่า อนุภาคจะอยู่บนพื้นผิวได้ก็ต่อเมื่ออนุภาคนั้นไม่ได้ถูกล้อมรอบด้วยอนุภาคอื่นๆ ทุกทิศทาง หรือ เราจะบอกว่ามันเป็นอนุภาคพื้นผิวถ้ามันไม่ถูกฝั่งลงไปในของกลุ่มก้อนนั้น เกณฑ์เหล่านี้ทำให้เราสามารถตัดสินใจว่ามันเป็นอนุภาคพื้นผิวหรือไม่ แต่เพื่อที่จะบอกคอมพิวเตอร์ให้เก็บข้อมูลของอนุภาคที่อยู่พื้นผิว เช่นตำแหน่งของอนุภาคเหล่านั้น เราก็ใช้หลักการที่จะอธิบายเหมือนกับหลักเกณฑ์ที่กล่าวมาก่อนหน้านี้

จุดหลายจุดที่เห็นในภาพด้านขวามือนั้นจะอยู่ในอนุภาค ซึ่งอนุภาคจะติดต่อกับอาณาเขตรูปกรวย ขอบเขตรูปกรวยนี้ถูกกำหนดขนาดพื้นที่ด้วยความยาว และมุมของมัน และเราเรียกอาณาเขตรูปกรวยที่ไม่มีอนุภาคใดๆ อยู่ภายในว่า "กรวยกลวง" (hollow cone) เรากำหนดให้อนุภาคอยู่บนพื้นผิวถ้าจุดอย่างน้อยหนึ่งจุดบนอนุภาคติดต่อกับกรวยกลวง
สำหรับกรวย เรามีตัวแปรที่เปลี่ยนแปลงได้ 2 ตัว คือ ความยาวของกรวย และขนาดของมุมกรวย ผลลัพธ์ที่แม่นยำหรือผลลัพธ์ที่น่าพอใจจะเกิดขึ้นได้ก็ขึ้นอยู่กับตัวแปร 2 ตัวนี้ ขั้นตอนวิธีแบบนี้สามารถปรับใช้ได้อย่างรวดเร็วกับสถานการณ์ที่มีโมเลกุลหลายโมเลกุลในโครงสร้าง หรือสถานการณ์ที่โมเลกุลนั้นมีรูอยู่ภายใน
== การทำงาน ==

ใช้วิธีการดัชนีเซลล์เพื่อเร่งให้เห็นอนุภาคข้างเคียง ความยาวด้านข้างของเซลล์กว้างกว่าความยาวด้านข้างของกรวยเล็กน้อย ซึ่งความยาวที่แตกต่างกันนี้จะถูกสังเกตและตัดออกไปในภายหลัง และอนุภาคข้างเคียงจะถูกระบุให้อยู่ภายในระยะที่ถูกตัดออกไปนั้น
ในการที่จะดูว่าอนุภาคที่เราสนใจมีกรวยกลวงหรือไม่ เราจะวนหาค่าผลรวมของอนุภาคข้างเคียง 3 อนุภาค เพื่อหาว่ามีอนุภาคข้างเคียงอื่นอยู่อีกหรือไม่ ผลรวมเชิงสเกล่าร์ของเวกเตอร์ที่เชื่อมต่อกับอนุภาคต่างๆ เหล่านี้ถูกใช้ในการพิจารณาว่าถ้าอาณาเขตรูปกรวยขยายไป 3 ช่วงอนุภาคแล้ว จะมีพื้นที่กว้างพอที่จะใส่ตัวแปร 2 ตัวหรือไม่ และมีอนุภาคอื่นอยู่ในอาณาเขตเขตรูปกรวยหรือไม่
ในการที่จะหาอนุภาคพื้นผิว ก็เพียงแค่หาอนุภาคที่มีกรวยกลวงอยู่ แม้จะมีกรวยกลวงเพียงอันเดียวก็จะถือว่าเป็นอนุภาคพื้นผิว และการวนเพื่อหาคำตอบจะสิ้นสุดลง
ในทางกลับกันเพื่อที่จะหาอนุภาคภายในพื้นผิวจำเป็นต้องค้นหาทุกๆ ค่าผลรวมของอนุภาคข้างเคียง 3 อนุภาค เพื่อที่จะยืนยันว่าไม่มีกรวยกลวงอยู่จริงๆ
แม้วิธีการง่ายๆ นี้จะเสร็จสมบูรณ์แต่เวลานั้นถูกใช้จนหมด ดังนั้นด้านล่างนี้จะกล่าวถึงวิธีอื่นๆ ที่ทำให้ความเร็วในการคำนวณของโปรแกรมดีขึ้นอย่างเห็นได้ชัด เพื่อจะลดเวลาการทำงานไม่ให้เวลาถูกใช้จนหมดไป

ในความเป็นจริง เราสามารถระบุหาตำแหน่งของอนุภาคภายในส่วนใหญ่ที่ห่างจากพื้นผิวได้อย่างง่ายดาย ตามนิยามที่ว่าพวกมันถูกล้อมรอบด้วยอนุภาคอื่นๆ จากรูปแสดงภาพจำลอง 2 มิติ เซลล์สีแดงเป็นอนุภาคภายใน เพราะถูกล้อมรอบด้วยเซลล์ข้างเคียง (สีน้ำเงิน) ถ้าเรามีจำนวนอนุภาคเยอะกว่านี้ต้องใช้แบบจำลอง 3 มิติในการพิจารณา
อนุภาคพื้นผิวบางอนุภาคสามารถหาได้ง่ายโดยวิธีทางเรขาคณิตที่ใช้คำนวณกลุ่มก้อนที่เป็นรูปทรงกลม ตัวอย่างเช่น ถ้าสมมุติให้กลุ่มก้อนเป็นรูปทรงกลม จะมีโอกาสที่จะเจอกรวยกลวง ได้มากขึ้น โดยการเลือกเวอร์เตอร์จากศูนย์กลางของกลุ่มก้อนถึงอนุภาคที่เราสนใจ ตามกฎทางเรขาคณิตของกรวย
=== การเกิดชั้นย่อยของขั้นตอนวิธีแบบกรวย ===

การเรียกใช้ขั้นตอนวิธีแบบกรวยหลายๆ ครั้งทำให้เกิดชั้นของพื้นผิวย่อยๆ ที่แตกต่างกันหลายชั้น โปรแกรมภาษาซีนั้นเป็นโปรแกรมที่ง่ายสำหรับการวิเคราะห์ชั้นพื้นผิวย่อยเหล่านี้ ซึ่งเกิดอยู่ในก้อนสีทองรูปหกเหลี่ยม (icosahedral gold cluster) ที่ประกอบด้วยอะตอม 2,624 อะตอม รูปนี้แสดงให้เห็นถึงชั้นย่อยของพื้นผิวในก้อนหกเหลี่ยม (แต่ละอะตอมในก้อนหกเหลี่ยมนี้เป็นอะตอมสีทองเหมือนกัน แต่ถูกระบายสีที่ต่างกันตามชั้นที่พวกมันสังกัดอยู่ )
== การเขียนโปรแกรมของขั้นตอนวิธีแบบกรวย ==
แม้ว่าภาษาซี (ซึ่งมีโครงสร้างภาษาแบบเชิงวัตถุ) จะเป็นภาษาที่เหมาะจะใช้เขียนขั้นตอนวิธีแบบนี้ แต่การโปรแกรมนั้นสามารถใช้กับโครงสร้างภาษาได้หลากหลายมากกว่าโครงสร้างแบบเชิงวัตถุอย่างเดียวการเขียนโปรแกรมควรจะสามารถใชักับข้อมูลขาเข้าและขาออกที่มีมาตรฐานตั้งแต่น้อยๆไปจนมีมาตรฐานมากๆได้ source code ของโปรแกรมควรจะอ่านง่าย และควรจะคอมเมนต์ที่เพียงพอเพื่อที่จะทำให้อ่านเข้าใจ
== รหัสในภาษาซี ==
=== การเรียกใช้ฟังก์ชัน ===
include 
include 
include 
include 
include "vector.h"
using namespace std;
=== การประกาศตัวแปรคลาส ===
double cosine; // cosine of the angle of the die cone
double cutoff2; // =cutoff*cutoff, cutoff is the side length of the die cone
double bin; // side length of cell = 1.0001 * cutoff
int num; // number of particles
Vector *q; // particle positions
int *d; // d[i]: 1--surface, -1--internal, 0--undetermined
int ***cell; // head particle of cell chain
int ***curr; // current-position pointer of cell chain
int nx, ny, nz; // cell numbers in 3 dimensions
int *list; // single chain of the particle list int the cells, end with -1
const double TINY = 1.0e-20; // zero value
void readConf;
void initCell;
void end;
void writeXYZ;
void markCell;
void markSurface;
int cmCheck( const int p, const int i, const int j, const int k );
int cmCone( const int p, const int np );
int cone( const Vector &v, const int p, const int np );
int realCheck( const int p, const int i, const int j, const int k );
bool realCone( const int p, const int np1, const int np2,
const int np3, const int *nbr, const int nnbr );
bool coneVector( const int p, const int np1, const int np2, const int np3,
Vector &v );
bool inLine( const int p, const int p1, const int p2 );
=== ฟังก์ชันต่างๆ ===
//---------------------------function readConf--------------------------
// Read one configuration of a molecule from the standard input.
// The input data should be in .xyz format.
void readConf{
cin >> num;
char s[1024];
cin.getline( s, 1024 );       // comment line
cin.getline( s, 1024 );
d =new int[num];
list = new int[num];
for( int i=0; i> s >> q[i];
} 
//------------------------------function initCell--------------------------------
// Find the boundary of the simulation box and locate particles into cells.
void initCell{
// Create cell array
Vector min( q[0].x, q[0].y, q[0].z );
Vector max( min );
for( int i=1; i max.x ) max.x = q[i].x;
if( q[i].y  max.y ) max.y = q[i].y;
if( q[i].z  max.z ) max.z = q[i].z;
}
nx = (int)ceil( (max.x-min.x)/bin );
ny = (int)ceil( (max.y-min.y)/bin );
nz = (int)ceil( (max.z-min.z)/bin );
cell = new int **[nx];
curr = new int **[nx];
for( int i=0; i=nx || y>=ny || z>=nz ||x
//------------------------- function end----------------------------
// Finishing job.
void end{
delete []q; delete []d;   delete []list;
for( int i=0; i
//-------------------------- function writeXYZ---------------------------
// Write out the configuration with the surface particles as Cd, the internal
// particles as Au. If any particles left undermined, they are maked as Ag.
void writeXYZ{
cout 
//---------------------------function markCell--------------------------
// Mark the particles in the internal cells. If a cell is internal is
// determined by looking around the designated neighbor cells to see if
// they all contain particles. A batch way to find the most interior particles.
void markCell{
for( int i=2; i
//---------------------------function markSurface--------------------------
// Mark the surface particles by looking at the bond angles with the neighbor
// particles and trying to find a hollow (without particles) cone region with
// angle larger than the designated threshold.
void markSurface{
// First round quick check for surface particles by center-of-mass normal
// vector criterion.
for( int i=0; i
//-------------------------function cmCheck----------------------------
// Check if a given particle is on the surface by cmCone.
int cmCheck( const int p, const int i, const int j, const int k ){
// Loop over all particles in this cell and the neighboring cells
for( int ni=i-1; ni=nx ) continue;
for( int nj=j-1; nj=ny ) continue;
for( int nk=k-1; nk=nz ) continue;
// internal particles are definitely out of the cone region
if( curr[ni][nj][nk] == -1 ) continue;
int np = cell[ni][nj][nk];
while( np != -1 ){
if( cmCone( p, np ) == 1 ) return 0;
np = list[ np ];
}
}
}
}
return 1;        // surface particle
}
//--------------------------- function cmCone--------------------------
// Check if np is in the cone region of p with the vector from center of mass
// to p and angle by cosine. If so, return true. Otherwise, return false.
// A quick way to pick up some surface particles before complete searching algorithm.
int cmCone( const int p, const int np ){
return cone( q[p].normalize, p, np );
}//-------------------------- function cone---------------------------
// Check if a given particle np is in the cone region centered at particle p
// with center vector of (nx,ny,nz) and cone angle theta
// with cos(theta)=cosine. Cone center vector must be normalized.
// Return 0 if the particle is not inside the region(s), 1 if in the cone
// region, -1 if in the inverse cone region.
int cone( const Vector &v, const int p, const int np ){
Vector d = (q[np] - q[p]).normalize;
double s = v*d;
if( s > cosine ) return 1;    // less than the angle
else if( s
//-------------------------- function realCheck---------------------------
// Check if the given particle is a surface particle without ambiguity by complete
// searching for a cone region larger than the die.
int realCheck( const int p, const int i, const int j, const int k ){
int *nbr = new int[1000];   // list of neighboring particles
int nnbr=0;    // number of neighbor particles
int *snbr = new int[1000];  // neighboring surface particles
int snnbr = 0;
int *unbr = new int[1000];  // neighboring undetermined particles
int unnbr = 0;
// Loop over all particles in this cell and the neighboring cells to make
// a list of neighboring particles.
for( int ni=i-1; ni=nx ) continue;
for( int nj=j-1; nj=ny ) continue;
for( int nk=k-1; nk=nz ) continue;
int np = cell[ni][nj][nk];
while( np != -1 ){
if( np != p )  {
// check if np is inside the cutoff region
Vector v = q[np] - q[p];
if( v.r2 
//-------------------------- function realCone---------------------------
// Real cone searching to find out if the cone region defined by
// (np1-p), (np2-p), (np3-p) or its inverse is hollow.
// Return true if at least one is hollow.
bool realCone( const int p, const int np1, const int np2,
const int np3, const int *nbr, const int nnbr ){
Vector v;
// Find the center vector of the cone spanned by ( np1, np2, np3 )
// centered at p.
if( !coneVector( p, np1, np2, np3, v ) ) return false;
// If any particles inside the cone region and the inverse cone region.
bool has = false;
bool hasInv = false;
for( int m=0; m
//--------------------------function coneVector---------------------------
// Find the center vector of the cone defined by the three vectors
// (np1-p), (np2-p), (np3-p).
// Return false if any error or the angle spanned is not big enough.
bool coneVector( const int p, const int np1, const int np2, const int np3,
Vector &v ){
const Vector v1 = (q[np1] - q[p]).normalize;
const Vector v2 = (q[np2] - q[p]).normalize;
const Vector v3 = (q[np3] - q[p]).normalize;
// it's the normal vector of plane spanned by v1, v2 and v3
v = (v1-v2)^(v2-v3);
// check if three points are on one line
double r = sqrt( v.r2 );
if( r  cosine ) // smaller angle
return false;
return true;
}
=== การเรียกใช้ฟังก์ชันหลัก ===
int main( int argc, char *argv[] ){
if( argc 
== อ้างอิง ==
Cone Algorithm By Yanting Wang
กขั้นตอนวิธีแบบกรวย
กขั้นตอนวิธีแบบกรวย