1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * http://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: See CVS logs. Details at http://www.graphviz.org/
9 *************************************************************************/
10
11#include "config.h"
12
13#ifdef HAVE_ANN
14//----------------------------------------------------------------------
15// File: get_nearest_neighb_graph.cpp
16//----------------------------------------------------------------------
17
18#include <cstdlib> // C standard library
19#include <cstdio> // C I/O (for sscanf)
20#include <cstring> // string manipulation
21#include <fstream> // file I/O
22#include <ANN/ANN.h> // ANN declarations
23
24using namespace std; // make std:: accessible
25int dim = 4; // dimension
26
27
28static void printPt(ostream &out, ANNpoint p) // print point
29{
30 out << "" << p[0];
31 for (int i = 1; i < dim; i++) {
32 out << "," << p[i];
33 }
34 out << "";
35}
36
37static void sortPtsX(int n, ANNpointArray pts){
38 /* sort so that edges always go from left to right in x-doordinate */
39 ANNpoint p;
40 ANNcoord x, y;
41 int i, j;
42 for (i = 0; i < n; i++){
43 for (j = 0; j < dim; j++){
44 p = pts[i];
45 if (p[0] < p[2] || (p[0] == p[2] && p[1] < p[3])) continue;
46 x = p[0]; y = p[1];
47 p[0] = p[2];
48 p[1] = p[3];
49 p[2] = x;
50 p[3] = y;
51 }
52 }
53}
54
55static void sortPtsY(int n, ANNpointArray pts){
56 /* sort so that edges always go from left to right in x-doordinate */
57 ANNpoint p;
58 ANNcoord x, y;
59 int i, j;
60 for (i = 0; i < n; i++){
61 for (j = 0; j < dim; j++){
62 p = pts[i];
63 if (p[1] < p[3] || (p[1] == p[3] && p[0] < p[2])) continue;
64 x = p[0]; y = p[1];
65 p[0] = p[2];
66 p[1] = p[3];
67 p[2] = x;
68 p[3] = y;
69 }
70 }
71}
72
73
74extern "C" void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0);
75
76void nearest_neighbor_graph_ann(int nPts, int dim, int k, double eps, double *x, int *nz0, int **irn0, int **jcn0, double **val0){
77
78 /* Gives a nearest neighbor graph is a list of dim-dimendional points. The connectivity is in irn/jcn, and the distance in val.
79
80 nPts: number of points
81 dim: dimension
82 k: number of neighbors needed
83 eps: error tolerance
84 x: nPts*dim vector. The i-th point is x[i*dim : i*dim + dim - 1]
85 nz: number of entries in the connectivity matrix irn/jcn/val
86 irn, jcn: the connectivity
87 val: the distance
88
89 note that there could be repeates
90 */
91
92 ANNpointArray dataPts; // data points
93 ANNidxArray nnIdx; // near neighbor indices
94 ANNdistArray dists; // near neighbor distances
95 ANNkd_tree* kdTree; // search structure
96
97 double *xx;
98 int *irn, *jcn;
99 double *val;
100 int nz;
101
102 irn = *irn0;
103 jcn = *jcn0;
104 val = *val0;
105
106
107 dataPts = annAllocPts(nPts, dim); // allocate data points
108 nnIdx = new ANNidx[k]; // allocate near neigh indices
109 dists = new ANNdist[k]; // allocate near neighbor dists
110
111 for (int i = 0; i < nPts; i++){
112 xx = dataPts[i];
113 for (int j = 0; j < dim; j++) xx[j] = x[i*dim + j];
114 }
115
116 //========= graph when sort based on x ========
117 nz = 0;
118 sortPtsX(nPts, dataPts);
119 kdTree = new ANNkd_tree( // build search structure
120 dataPts, // the data points
121 nPts, // number of points
122 dim); // dimension of space
123 for (int ip = 0; ip < nPts; ip++){
124 kdTree->annkSearch( // search
125 dataPts[ip], // query point
126 k, // number of near neighbors
127 nnIdx, // nearest neighbors (returned)
128 dists, // distance (returned)
129 eps); // error bound
130
131 for (int i = 0; i < k; i++) { // print summary
132 if (nnIdx[i] == ip) continue;
133 val[nz] = dists[i];
134 irn[nz] = ip;
135 jcn[nz++] = nnIdx[i];
136 //cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n";
137 //printPt(cout, dataPts[ip]);
138 //cout << "--";
139 //printPt(cout, dataPts[nnIdx[i]]);
140 //cout << "\n";
141 }
142 }
143
144
145 //========= graph when sort based on y ========
146 sortPtsY(nPts, dataPts);
147 kdTree = new ANNkd_tree( // build search structure
148 dataPts, // the data points
149 nPts, // number of points
150 dim); // dimension of space
151 for (int ip = 0; ip < nPts; ip++){
152 kdTree->annkSearch( // search
153 dataPts[ip], // query point
154 k, // number of near neighbors
155 nnIdx, // nearest neighbors (returned)
156 dists, // distance (returned)
157 eps); // error bound
158
159 for (int i = 0; i < k; i++) { // print summary
160 if (nnIdx[i] == ip) continue;
161 val[nz] = dists[i];
162 irn[nz] = ip;
163 jcn[nz++] = nnIdx[i];
164 // cout << ip << "--" << nnIdx[i] << " [len = " << dists[i]<< ", weight = \"" << 1./(dists[i]) << "\", wt = \"" << 1./(dists[i]) << "\"]\n";
165 }
166 }
167
168 delete [] nnIdx; // clean things up
169 delete [] dists;
170 delete kdTree;
171
172 *nz0 = nz;
173
174 annClose(); // done with ANN
175
176
177
178}
179
180#endif /* HAVE_ANN */
181