Index: Damkjer/Util/Math/par_eig.cpp =================================================================== --- Damkjer/Util/Math/par_eig.cpp (revision 8) +++ Damkjer/Util/Math/par_eig.cpp (revision 8) @@ -0,0 +1,154 @@ +//========================================================================= +// FILE: par_eig.cpp +// +// Copyright (C) 2012 Kristian Damkjer. +// +// DESCRIPTION: This MEX source file provides a fast implementation of cov +// method for cell-arrays of real valued matrices. +// +// LIMITATIONS: Does not work for cell-arrays of complex matrices. +// +// SOFTWARE HISTORY: +//> 2013-JUL-03 K. Damkjer +// Initial Coding. +//< +//========================================================================= + +#ifdef _OPENMP +#include +#endif + +#include +#include + +#if MATLAB_MAJOR <= 7 && MATLAB_MINOR <= 10 && defined(_CHAR16T) + #define CHAR16_T +#endif + +#include + +#include "mex.h" + +void mexFunction( + int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + if (nrhs != 1 || !mxIsCell(prhs[0])) + { + mexErrMsgIdAndTxt("Damkjer:fastcov:varargin", + "Missing or invalid input argument."); + } + + if (nlhs > 2) + { + mexErrMsgIdAndTxt("Damkjer:fastcov:varargout", + "Too many output arguments."); + } + + mwSize cells = mxGetNumberOfElements (prhs[0]); + + plhs[0] = mxCreateCellMatrix(cells, 1); + + // Better way? + if (nlhs > 1) + { + plhs[1] = mxCreateCellMatrix(cells, 1); + } + + std::vector data(cells,0); + std::vector Ms(cells,0); + std::vector Ns(cells,0); + + std::vector vecs(cells,0); + std::vector vecs_data(cells,0); + + std::vector vals(cells,0); + std::vector vals_data(cells,0); + + // Note for future: Ms - points, Ns - dimensions + for (int cell = 0; cell < cells; ++cell) + { + data[cell]=mxGetPr(mxGetCell(prhs[0], cell)); + Ms[cell]=mxGetM(mxGetCell(prhs[0], cell)); + Ns[cell]=mxGetN(mxGetCell(prhs[0], cell)); + } + + // Always return values as vector + for (int cell = 0; cell < cells; ++cell) + { + // We will be setting each value, so don't bother to initialize to zero. + vals[cell] = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(vals[cell], Ms[cell]); + mxSetN(vals[cell], 1); + mxSetData(vals[cell], mxMalloc(sizeof(double)*Ms[cell])); + vals_data[cell] = mxGetPr(vals[cell]); + } + +#ifdef _OPENMP + omp_set_dynamic(1); + omp_set_num_threads(omp_get_num_procs()); +#endif + + #pragma omp parallel for + for (int cellp = 0; cellp < cells; ++cellp) + { + Eigen::VectorXd eivals = + Eigen::Map(data[cellp], Ms[cellp], Ns[cellp]). + selfadjointView().eigenvalues(); + + for (mwSize m = Ms[cellp]; m --> 0;) + { + vals_data[cellp][m] = eivals(m); + } + } + + if (nlhs < 2) + { + // Only return values as vector + for (int cell = 0; cell < cells; ++cell) + { + mxSetCell(plhs[0], cell, vals[cell]); + } + } + else + { + // Return both values and vectors as matices + for (int cell = 0; cell < cells; ++cell) + { + // We will be setting each value, so don't bother to initialize to zero. + vecs[cell] = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(vecs[cell], Ms[cell]); + mxSetN(vecs[cell], Ns[cell]); + mxSetData(vecs[cell], mxMalloc(sizeof(double)*Ms[cell]*Ns[cell])); + vecs_data[cell] = mxGetPr(vecs[cell]); + } + +#ifdef _OPENMP + omp_set_dynamic(1); + omp_set_num_threads(omp_get_num_procs()); +#endif + + #pragma omp parallel for + for (int cellp = 0; cellp < cells; ++cellp) + { + Eigen::SelfAdjointEigenSolver + eig(Eigen::Map(data[cellp], Ms[cellp], Ns[cellp])); + + Eigen::MatrixXd eivecs=eig.eigenvectors(); + + for (mwSize m = Ms[cellp]; m --> 0;) + { + for (mwSize n = Ns[cellp]; n --> 0;) + { + vecs_data[cellp][n + Ns[cellp] * m] = eivecs(n,m); + } + } + } + + for (int cell = 0; cell < cells; ++cell) + { + mxSetCell(plhs[0], cell, vecs[cell]); + mxSetCell(plhs[1], cell, vals[cell]); + } + } +} Index: Damkjer/Util/SpatialIndexing/VpTree/VpTreeAPI.cpp =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/VpTreeAPI.cpp (revision 8) +++ Damkjer/Util/SpatialIndexing/VpTree/VpTreeAPI.cpp (revision 8) @@ -0,0 +1,752 @@ +//***************************************************************************** +// FILE: VpTreeAPI.cpp +// +// Copyright (C) 2012 Kristian Damkjer. +// +// DESCRIPTION: This MATLAB Executable (MEX) gateway provides the application +// interface for working with VP Trees in MATLAB. All functional +// interfaces are routed through this single function to ensure +// locks and memory management are handled properly by MATLAB. +// +// LIMITATIONS: All interfaces to the VP Tree objects must be routed through +// this gateway file. This interface is required until and unless +// MATLAB provides a mechanism for modifying locks on MEX files +// other than the current file. +// +// SOFTWARE HISTORY: +//> 2012-SEP-11 K. Damkjer +// Initial Coding. +// 2013-JUL-19 K. Damkjer +// Improved structure and documentation. Prevented double free of +// VP Tree memory and infinite locking of MEX library. +//< +//***************************************************************************** + +#ifdef _OPENMP +#include +#endif + +#include "Util/MATLAB/ClassHandle.h" +#include "VpTree.h" + +// Fix "wide char" definition for older versions of MATLAB. This must be placed +// after other includes and before the mex.h include. +#if MATLAB_MAJOR <= 7 && MATLAB_MINOR <= 10 && defined(_CHAR16T) + #define CHAR16_T +#endif + +#include "mex.h" + +namespace Damkjer +{ +//***************************************************************************** +// FUNCTION: vpTreeCreate +//***************************************************************************** +void vpTreeCreate(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]); +//> This function provides the interface to the VpTree contructors. New objects +// are wrapped by the ClassHandle template to ensure that locks are managed +// correctly to avoid double frees of memory managed by this class, to ensure +// objects have a life span separate from this function's scope, and to +// provide convenience methods for passing the pointer to the object between +// MATLAB and C++. +// +// See the VpTree class documentation for details on the method. +// +// The function signature is identical to the mexFunction gateway for +// convenience only. +//< + +//***************************************************************************** +// FUNCTION: vpTreeDestroy +//***************************************************************************** +void vpTreeDestroy(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]); +//> This function provides the interface to the VpTree ClassHandle wrapper +// destructor. Destroying objects through the ClassHandle wrapper ensures +// memory leaks are not introduced through this MEX function and allows the +// compiled MEX library to be unlocked once all objects have been destroyed +// and their memory returned to the system. +// +// See the VpTree class documentation for details on the method. +// +// The function signature is identical to the mexFunction gateway for +// convenience only. +//< + +//***************************************************************************** +// FUNCTION: vpTreeFRANN +//***************************************************************************** +void vpTreeFRANN(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]); +//> This function provides the interface to the VpTree fixed-radius all nearest +// neighbor (FRANN) search. +// +// See the VpTree class documentation for details on the method. +// +// The function signature is identical to the mexFunction gateway for +// convenience only. +//< + +//***************************************************************************** +// FUNCTION: vpTreeKANN +//***************************************************************************** +void vpTreeKANN(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]); +//> This function provides the interface to the VpTree k all nearest neighbor +// (KANN) search. +// +// See the VpTree class documentation for details on the method. +// +// The function signature is identical to the mexFunction gateway for +// convenience only. +//< +} + +//***************************************************************************** +// FUNCTION: mexFunction +// +//> The MATLAB Executable Gateway Function. +//< +//***************************************************************************** +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //*** + // This MEX function should always have at least one input parameter to + // handle function dispatching through the MEX gateway. + //*** + if (nrhs < 1) + { + mexErrMsgIdAndTxt("Damkjer:VpTreeAPI:nargin", + "VpTreeAPI requires at least one input which " + "specifies the operation to be performed."); + } + + // The operation switch is simply a string. + char* operation = mxArrayToString(prhs[0]); + + if (!operation) + { + mexErrMsgIdAndTxt("Damkjer:VpTreeAPI:InvalidOperation", + "Invalid mode supplied to VpTreeAPI."); + } + + // Dispatch to helper functions. Err if operation is not recognized. + if (!strcmp("create", operation)) + { + Damkjer::vpTreeCreate(nlhs, plhs, nrhs, prhs); + } + else if (!strcmp("destroy", operation)) + { + Damkjer::vpTreeDestroy(nlhs, plhs, nrhs, prhs); + } + else if (!strcmp("search_frann", operation)) + { + Damkjer::vpTreeFRANN(nlhs, plhs, nrhs, prhs); + } + else if (!strcmp("search_kann", operation)) + { + Damkjer::vpTreeKANN(nlhs, plhs, nrhs, prhs); + } + else + { + mexErrMsgIdAndTxt("Damkjer:VpTreeAPI:UnknownOperation", + "Unrecognized mode provided to VpTreeAPI."); + } + + // Prevent a slow memory leak. + mxFree(operation); +} + +namespace Damkjer +{ +//***************************************************************************** +// FUNCTION: vpTreeCreate +//***************************************************************************** +void vpTreeCreate(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs < 2 || nrhs > 2) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeCreate:nargin", + "The vpTreeCreate function requires a single input."); + } + + if (nlhs > 1) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeCreate:nargout", + "The vpTreeCreate function requires a single output."); + } + + //*** + // The single input required by vpTreeCreate is the collection of points to + // be indexed. + //*** + const mxArray* points=prhs[1]; + + // Check to make sure that points are real-valued numerics + if (mxIsSparse(points) || + mxGetNumberOfDimensions(points) != 2 || + mxIsComplex(points) || + !mxIsNumeric(points)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeCreate:prhs", + "Input to vpTreeCreate must be a full, 2-D matrix " + "of real-valued data representing N-dimensional " + "observations."); + } + + // Attempt to make the VP Tree. + const mwSize dims = mxGetM(points); + const mwSize elems = mxGetN(points); + + double* data = mxGetPr(points); + + //*** + // Selection of data structure for the point database is not arbitrary. + // + // 1. Point coordinates are passed in as double-precision values. Preserve + // the data fidelity. + // 2. Individual points are often low-dimensional. Keep the coordinates for + // a single point in contiguous memory by using a std::vector container. + // 3. The point database can be quite large. Since VP Tree is an online data + // structure, allow references to points to occupy non-contiguous memory. + // This also speeds construction and destruction of the point database. + //*** + std::deque > pointData(elems, + std::vector(dims)); + + for (mwSize elem = elems; elem --> 0;) + { + for (mwSize dim = dims; dim --> 0;) + { + pointData[elem][dim]=data[elem*dims+dim]; + } + } + + //*** + // Use the ClassHandle named constructor to provide a handle to the new VP + // Tree object. + //*** + plhs[0] = ptrAsMat(new VpTree<>(pointData)); +} + +//***************************************************************************** +// FUNCTION: vpTreeDestroy +//***************************************************************************** +void vpTreeDestroy(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs != 2 || !mxIsNumeric(prhs[1])) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeDestroy:varargin", + "The vpTreeDestroy function requires a single " + "input."); + } + + // Free the VP Tree through the ClassHandle named destructor. + destroyHandleTo >(prhs[1]); +} + +//***************************************************************************** +// FUNCTION: vpTreeFRANN +//***************************************************************************** +void vpTreeFRANN(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs < 3 || nrhs > 4 || !mxIsNumeric(prhs[1])) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeFRANN:varargin", + "Invalid number of input arguments."); + } + + if (nlhs>2) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeFRANN:varargout", + "Invalid number of output arguments."); + } + + // Retrieve the tree object through the ClassHandle helper method. + const VpTree<>& tree = matAsObj >(prhs[1]); + + // The second parameter should be the query points. + const mxArray* queries=prhs[2]; + + // Check to make sure that query points are real-valued numerics + if (mxIsSparse(queries) || + mxGetNumberOfDimensions(queries) != 2 || + mxIsComplex(queries) || + !mxIsNumeric(queries)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeFRANN:prhs", + "Second input parameter to vpTreeFRANN must be a " + "full, 2-D matrix of real-valued data representing " + "multi-dimensional queries."); + } + + // The third parameter should be the query radius. + const mxArray* rData=prhs[3]; + + // Check to make sure that radius is a real-valued numeric scalar. + if (mxIsSparse(rData) || + mxGetNumberOfElements(rData) != 1 || + mxIsComplex(rData) || + !mxIsNumeric(rData)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeFRANN:prhs", + "Third input parameter to vpTreeFRANN must be a " + "real-valued scalar representing desired neighborhood " + "radius limit."); + } + + //*** + // Get the query points. + // + // Selection of data structure for the query database is not arbitrary. + // + // 1. Point coordinates are passed in as double-precision values. Preserve + // the data fidelity. + // 2. Individual points are often low-dimensional. Keep the coordinates for + // a single point in contiguous memory by using a std::vector container. + // 3. The query database can be quite large. Allow references to points to + // occupy non-contiguous memory. + //*** + const mwSize dims = mxGetM(queries); + const mwSize elems = mxGetN(queries); + + double* data = mxGetPr(queries); + std::deque > queryData(elems, std::vector(dims)); + + for (mwIndex elem = elems; elem --> 0;) + { + for (mwIndex dim = dims; dim --> 0;) + { + queryData[elem][dim]=data[elem*dims+dim]; + } + } + + // Get the desired neighborhood radius limit. + double radius = (*(double*)mxGetData(rData)); + + //*** + // The first output parameter holds the nearest neighbor indices. This + // collection is represented as a cell array of vectors with a cell for each + // query point. + //*** + plhs[0] = mxCreateCellMatrix(elems, 1); + + //*** + // If desired, the distance from the query point to each of the nearest + // neighbors can also be provided. This collection is represented as a cell + // array of vectors with a cell for each query point. + //*** + if (nlhs==2) + { + plhs[1] = mxCreateCellMatrix(elems, 1); + } + + //*** + // The following logic flows may seem counter-intuitive, but are structured + // to intentionally avoid nested critical sections in parallelized code. + // + // Unfortunately, any calls to the MATLAB API must be treated as belonging + // to a critical section since none of the API is thread-safe. + //*** + +#ifdef _OPENMP + omp_set_dynamic(1); + omp_set_num_threads(omp_get_num_procs()); +#endif + + //*** + // The VP Tree data structure always provides results as a set of pairs of + // indices and distances. There is no performance benefit to omitting the + // distance computations since they are required for the search algorithm. + //*** + std::vector, + std::deque > > results(queryData.size()); + + //*** + // The first embarassingly parallelizable section simply searches for each + // query point's neighbors in parallel. This is much simpler than attempting + // to parallelize the search for a single point, and likely to yield + // superior results. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + results[q] = tree.rnn(queryData[q], radius); + } + + //*** + // Allocating memory for the results to be passed back to MATLAB must take + // place in a "critical section" since it involves exercising the MEX API. + //*** + std::vector nbr_idxs; + std::vector point_idxs; + + nbr_idxs.reserve(queryData.size()); + point_idxs.reserve(queryData.size()); + + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + nbr_idxs.push_back(mxCreateNumericMatrix(0, 0, mxINDEX_CLASS, mxREAL)); + mxSetM(nbr_idxs[q], neighbors); + mxSetN(nbr_idxs[q], 1); + mxSetData(nbr_idxs[q], mxMalloc(sizeof(mwIndex)*neighbors*1)); + + point_idxs.push_back((mwIndex*)mxGetData(nbr_idxs[q])); + } + + //*** + // Once memory has been allocated, the actual results can be populated in + // parallel. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + mxArray* neighbor_idxs = nbr_idxs[q]; + mwIndex* idxs = point_idxs[q]; + + for (mwIndex idx = neighbors; idx --> 0;) + { + idxs[idx]=results[q].first[idx]+1; + } + } + + //*** + // Marking the data for return to MATLAB must take place in a "critical + // section" since it involves exercising the MEX API. This also changes + // ownership and memory management responsibilities to MATLAB. We will not + // free this data. + //*** + for (int q = 0; q < queryData.size(); ++q) + { + mxSetCell(plhs[0], q, nbr_idxs[q]); + } + + // Repeat the "hand-off" to MATLAB for distance data, if it was requested. + if (nlhs==2) + { + //*** + // Allocating memory for the results to be passed back to MATLAB must + // take place in a "critical section" since it involves exercising the + // MEX API. + //*** + std::vector nbr_dists; + std::vector point_dists; + + nbr_dists.reserve(queryData.size()); + point_dists.reserve(queryData.size()); + + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + nbr_dists[q] = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(nbr_dists[q], neighbors); + mxSetN(nbr_dists[q], 1); + mxSetData(nbr_dists[q], mxMalloc(sizeof(double)*neighbors*1)); + + point_dists.push_back(mxGetPr(nbr_dists[q])); + } + + //*** + // Once memory has been allocated, the actual results can be populated in + // parallel. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + mxArray* neighbor_dists = nbr_dists[q]; + double* dists = point_dists[q]; + + for (mwIndex idx = neighbors; idx --> 0;) + { + dists[idx]=results[q].second[idx]; + } + } + + //*** + // Marking the data for return to MATLAB must take place in a "critical + // section" since it involves exercising the MEX API. This also changes + // ownership and memory management responsibilities to MATLAB. We will + // not free this data. + //*** + for (int q = 0; q < queryData.size(); ++q) + { + mxSetCell(plhs[1], q, nbr_dists[q]); + } + } +} + +//***************************************************************************** +// FUNCTION: vpTreeKANN +//***************************************************************************** +void vpTreeKANN(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs != 5 || !mxIsNumeric(prhs[1])) + { + mexErrMsgIdAndTxt("Damkjer:kannVpTree:varargin", + "Invalid number of arguments"); + } + + // Retrieve the tree object through the ClassHandle helper method. + const VpTree<>& tree = matAsObj >(prhs[1]); + + // The second parameter should be the query points. + const mxArray* queries=prhs[2]; + + // Check to make sure that query points are real-valued numerics. + if (mxIsSparse(queries) || + mxGetNumberOfDimensions(queries) != 2 || + mxIsComplex(queries) || + !mxIsNumeric(queries)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeKANN:prhs", + "Second input parameter to vpTreeKANN must be a full, " + "2-D matrix of real-valued data representing " + "multi-dimensional queries."); + } + + // The third parameter should be the desired neighborhood cardinality. + const mxArray* kData=prhs[3]; + + // Check to make sure that cardinality is a real-valued numeric scalar. + if (mxIsSparse(kData) || + mxGetNumberOfElements(kData) != 1 || + mxIsComplex(kData) || + !mxIsNumeric(kData)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeKANN:prhs", + "Third input parameter to vpTreeKANN must be an " + "real-valued scalar representing desired neighborhood " + "cardinality."); + } + + // The fourth parameter should be the desired radius limit. + const mxArray* rData=prhs[4]; + + // Check to make sure that radius is a real-valued numeric scalar. + if (mxIsSparse(rData) || + mxGetNumberOfElements(rData) != 1 || + mxIsComplex(rData) || + !mxIsNumeric(rData)) + { + mexErrMsgIdAndTxt("Damkjer:vpTreeKANN:prhs", + "Fourth input parameter to vpTreeKANN must be an " + "real-valued scalar representing desired neighborhood " + "radius limit."); + } + + //*** + // Get the query points. + // + // Selection of data structure for the query database is not arbitrary. + // + // 1. Point coordinates are passed in as double-precision values. Preserve + // the data fidelity. + // 2. Individual points are often low-dimensional. Keep the coordinates for + // a single point in contiguous memory by using a std::vector container. + // 3. The query database can be quite large. Allow references to points to + // occupy non-contiguous memory. + //*** + const mwSize dims = mxGetM(queries); + const mwSize elems = mxGetN(queries); + + double* data = mxGetPr(queries); + std::deque > queryData(elems, + std::vector(dims)); + + for (mwIndex elem = elems; elem --> 0;) + { + for (mwIndex dim = dims; dim --> 0;) + { + queryData[elem][dim]=data[elem*dims+dim]; + } + } + + // Get the desired neighborhood cardinality. + mwSize k = (mwSize)(*((double*)mxGetData(kData))); + + // Get the desired neighborhood radius limit. + double radius = (*(double*)mxGetData(rData)); + + //*** + // The first output parameter holds the nearest neighbor indices. This + // collection is represented as a cell array of vectors with a cell for each + // query point. + //*** + plhs[0] = mxCreateCellMatrix(elems, 1); + + //*** + // If desired, the distance from the query point to each of the nearest + // neighbors can also be provided. This collection is represented as a cell + // array of vectors with a cell for each query point. + //*** + if (nlhs==2) + { + plhs[1] = mxCreateCellMatrix(elems, 1); + } + + //*** + // The following logic flows may seem counter-intuitive, but are structured + // to intentionally avoid nested critical sections in parallelized code. + // + // Unfortunately, any calls to the MATLAB API must be treated as belonging + // to a critical section since none of the API is thread-safe. + //*** + +#ifdef _OPENMP + omp_set_dynamic(1); + omp_set_num_threads(omp_get_num_procs()); +#endif + + //*** + // The VP Tree data structure always provides results as a set of pairs of + // indices and distances. There is no performance benefit to omitting the + // distance computations since they are required for the search algorithm. + //*** + std::vector, + std::deque > > results(queryData.size()); + + //*** + // The first embarassingly parallelizable section simply searches for each + // query point's neighbors in parallel. This is much simpler than attempting + // to parallelize the search for a single point, and likely to yield + // superior results. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + results[q] = tree.knn(queryData[q], k, radius); + } + + //*** + // Allocating memory for the results to be passed back to MATLAB must take + // place in a "critical section" since it involves exercising the MEX API. + //*** + std::vector nbr_idxs; + std::vector point_idxs; + + nbr_idxs.reserve(queryData.size()); + point_idxs.reserve(queryData.size()); + + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + nbr_idxs.push_back(mxCreateNumericMatrix(0, 0, mxINDEX_CLASS, mxREAL)); + mxSetM(nbr_idxs[q], neighbors); + mxSetN(nbr_idxs[q], 1); + mxSetData(nbr_idxs[q], mxMalloc(sizeof(mwIndex)*neighbors*1)); + + point_idxs.push_back((mwIndex*)mxGetData(nbr_idxs[q])); + } + + //*** + // Once memory has been allocated, the actual results can be populated in + // parallel. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + mxArray* neighbor_idxs = nbr_idxs[q]; + mwIndex* idxs = point_idxs[q]; + + for (mwIndex idx = neighbors; idx --> 0;) + { + idxs[idx]=results[q].first[idx]+1; + } + } + + //*** + // Marking the data for return to MATLAB must take place in a "critical + // section" since it involves exercising the MEX API. This also changes + // ownership and memory management responsibilities to MATLAB. We will not + // free this data. + //*** + for (int q = 0; q < queryData.size(); ++q) + { + mxSetCell(plhs[0], q, nbr_idxs[q]); + } + + // Repeat the "hand-off" to MATLAB for distance data, if it was requested. + if (nlhs==2) + { + //*** + // Allocating memory for the results to be passed back to MATLAB must + // take place in a "critical section" since it involves exercising the + // MEX API. + //*** + std::vector nbr_dists; + std::vector point_dists; + + nbr_dists.reserve(queryData.size()); + point_dists.reserve(queryData.size()); + + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + nbr_dists[q] = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(nbr_dists[q], neighbors); + mxSetN(nbr_dists[q], 1); + mxSetData(nbr_dists[q], mxMalloc(sizeof(double)*neighbors*1)); + + point_dists.push_back(mxGetPr(nbr_dists[q])); + } + + //*** + // Once memory has been allocated, the actual results can be populated in + // parallel. + //*** + #pragma omp parallel for + for (int q = 0; q < queryData.size(); ++q) + { + mwSize neighbors = results[q].first.size(); + + mxArray* neighbor_dists = nbr_dists[q]; + double* dists = point_dists[q]; + + for (mwIndex idx = neighbors; idx --> 0;) + { + dists[idx]=results[q].second[idx]; + } + } + + //*** + // Marking the data for return to MATLAB must take place in a "critical + // section" since it involves exercising the MEX API. This also changes + // ownership and memory management responsibilities to MATLAB. We will + // not free this data. + //*** + for (int q = 0; q < queryData.size(); ++q) + { + mxSetCell(plhs[1], q, nbr_dists[q]); + } + } +} + +} Index: mkjer/Util/SpatialIndexing/VpTree/deleteVpTree.cpp =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/deleteVpTree.cpp (revision 0) +++ (revision ) @@ -1,41 +1,0 @@ -//************************************************************************* -// FILE: deleteVpTree.cpp -// -// Copyright (C) 2012 Kristian Damkjer. -// -// DESCRIPTION: This MEX source file safely deletes the object pointed to -// by the encoded raw pointer. -// -// LIMITATIONS: Unfortunately, the raw pointer trick is necessary to -// persist objects between MEX function calls. Ideally, we -// would simply return the created object. However, that would -// violate MATLAB's current internal memory management system. -// -// SOFTWARE HISTORY: -//> 2012-SEP-11 K. Damkjer -// Initial Coding. -//< -//************************************************************************* - -#include "Util/MATLAB/ClassHandle.h" -#include "VpTree.h" - -#ifdef _CHAR16T -#define CHAR16_T -#endif - -#include "mex.h" - -void mexFunction( - int nlhs, mxArray * plhs[], - int nrhs, const mxArray * prhs[]) -{ - // check the arguments - if (nrhs!=1 || !mxIsNumeric(prhs[0])) - { - mexErrMsgIdAndTxt("Damkjer:deleteVpTree:varargin", - "varargin{1} must be a valid VpTree pointer"); - } - - destroyHandleTo >(prhs[0]); -} Index: mkjer/Util/SpatialIndexing/VpTree/frannVpTree.cpp =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/frannVpTree.cpp (revision 0) +++ (revision ) @@ -1,173 +1,0 @@ -//========================================================================= -// FILE: frannVpTree.cpp -// -// Copyright (C) 2012 Kristian Damkjer. -// -// DESCRIPTION: This MEX source file provides an implementation of the -// fixed-radius all nearest neighbors search algorithm given -// a set of query points, a radius, and a VpTree representing -// the point data base. -// -// LIMITATIONS: This function has the potential to consume a large amount -// of memory. While I have confirmed that there are no memory -// leaks in the code to the best of my current ability, I have -// induced predictable crashes in MATLAB by finding all -// nearest neighbors for a 4M point data base. The crash did -// not come from this code, but rather from Handle Graphics -// when attempting to render a waitbar. As a result, it is -// recommended that points be streamed or blocked together -// into more manageable chunks. -// -// SOFTWARE HISTORY: -//> 2012-SEP-11 K. Damkjer -// Initial Coding. -//< -//========================================================================= - -#ifdef _OPENMP -#include -#endif - -#include "Util/MATLAB/ClassHandle.h" -#include "VpTree.h" - -#ifdef _CHAR16T -#define CHAR16_T -#endif - -#include "mex.h" - -void mexFunction( - int nlhs, mxArray* plhs[], - int nrhs, const mxArray* prhs[]) -{ - // check the arguments - if (nrhs<2 || nrhs > 3 || !mxIsNumeric(prhs[0])) - { - mexErrMsgIdAndTxt("Damkjer:deleteVpTree:varargin", - "varargin{1} must be a valid VpTree pointer"); - } - - if (nlhs>2) - { - mexErrMsgIdAndTxt("Damkjer:deleteVpTree:varargout", - "Invalid number of output arguments."); - } - - // retrieve the tree pointer - const VpTree<>& tree = matAsObj >(prhs[0]); - - const mxArray* points=prhs[1]; - - // Check to make sure that points are real numerics - if (mxIsSparse(points) || - mxGetNumberOfDimensions(points) != 2 || - mxIsComplex(points) || - !mxIsNumeric(points)) - { - mexErrMsgIdAndTxt("Damkjer:makeVpTree:prhs", - "Point input to kannVpTree must be a full, 2-D matrix representing ND queries."); - } - - const mwSize dims = mxGetM(points); - const mwSize elems = mxGetN(points); - - double* data = mxGetPr(points); - std::deque > pointData(elems, std::vector(dims)); - - for (mwIndex elem = elems; elem --> 0;) - { - for (mwIndex dim = dims; dim --> 0;) - { - pointData[elem][dim]=data[elem*dims+dim]; - } - } - - double radius; - - if (nrhs == 2) - { - radius = 1; - } - else - { - const mxArray* rData=prhs[2]; - - if (mxIsSparse(rData) || - mxGetNumberOfElements(rData) != 1 || - mxIsComplex(rData) || - !mxIsNumeric(rData)) - { - mexErrMsgIdAndTxt("Damkjer:makeVpTree:prhs", - "Input to makeVpTree must be a full, 2-D matrix representing ND observations."); - } - - radius = mxGetScalar(rData); - } - - plhs[0] = mxCreateCellMatrix(elems, 1); - - if (nlhs==2) - { - plhs[1] = mxCreateCellMatrix(elems, 1); - } - - #ifdef _OPENMP - omp_set_dynamic(1); - omp_set_num_threads(omp_get_num_procs()); - #endif - #pragma omp parallel for - for (int p = 0; p < pointData.size(); ++p) - { - std::vector q=pointData[p]; - - std::pair, - std::deque > results = tree.rnn(q, radius); - - mwSize neighbors = results.first.size(); - - mxArray* neighbor_idxs = 0; - mwIndex* idxs = 0; - mxArray* neighbor_dists = 0; - double* dists = 0; - - #pragma omp critical - { - neighbor_idxs = mxCreateNumericMatrix(neighbors, 1, mxINDEX_CLASS, mxREAL); - - idxs = (mwIndex*)mxGetData(neighbor_idxs); - } - - for (mwIndex idx = neighbors; idx --> 0;) - { - idxs[idx]=results.first[idx]+1; - } - - #pragma omp critical - { - mxSetCell(plhs[0], p, neighbor_idxs); - } - - if (nlhs==2) - { - #pragma omp critical - { - neighbor_dists = mxCreateDoubleMatrix(neighbors, 1, mxREAL); - - dists = mxGetPr(neighbor_dists); - } - - for (mwIndex idx = neighbors; idx --> 0;) - { - dists[idx]=results.second[idx]; - } - - #pragma omp critical - { - mxSetCell(plhs[1], p, neighbor_dists); - } - } - } - - return; -} Index: mkjer/Util/SpatialIndexing/VpTree/kannVpTree.cpp =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/kannVpTree.cpp (revision 5) +++ (revision ) @@ -1,193 +1,0 @@ -//========================================================================= -// FILE: kannVpTree.cpp -// -// Copyright (C) 2012 Kristian Damkjer. -// -// DESCRIPTION: This MEX source file provides an implementation of the -// k all nearest neighbors search algorithm given a set of -// query points, a number of neighbors, and a VpTree -// representing the point data base. -// -// LIMITATIONS: This function has the potential to consume a large amount -// of memory. While I have confirmed that there are no memory -// leaks in the code to the best of my current ability, I have -// induced predictable crashes in MATLAB in the related -// MEX function, frannVpTree, by finding all nearest neighbors -// for a 4M point data base. The crash did not come from this -// code, but rather from Handle Graphics when attempting to -// render a waitbar. As a result, it is recommended that -// points be streamed or blocked together into more manageable -// chunks. -// -// SOFTWARE HISTORY: -//> 2012-SEP-11 K. Damkjer -// Initial Coding. -// 2013-FEB-04 K. Damkjer -// Restructure parallel sections to see if any performance -// gain is realized. -//< -//========================================================================= - -#include - -#ifdef _OPENMP -#include -#endif - -#include "Util/MATLAB/ClassHandle.h" -#include "VpTree.h" - -#ifdef _CHAR16T -#define CHAR16_T -#endif - -#include "mex.h" - -void mexFunction( - int nlhs, mxArray* plhs[], - int nrhs, const mxArray* prhs[]) -{ - // check the arguments - if (nrhs != 4 || !mxIsNumeric(prhs[0])) - { - mexErrMsgIdAndTxt("Damkjer:kannVpTree:varargin", - "Invalid number of arguments"); - } - - // retrieve the tree pointer - const VpTree<>& tree = matAsObj >(prhs[0]); - - const mxArray* points=prhs[1]; - - // Check to make sure that points are real numerics - if (mxIsSparse(points) || - mxGetNumberOfDimensions(points) != 2 || - mxIsComplex(points) || - !mxIsNumeric(points)) - { - mexErrMsgIdAndTxt("Damkjer:kannVpTree:prhs", - "Point input to kannVpTree must be a full, 2-D matrix representing ND queries."); - } - - const mwSize dims = mxGetM(points); - const mwSize elems = mxGetN(points); - - double* data = mxGetPr(points); - std::deque > pointData(elems, std::vector(dims)); - - for (mwIndex elem = elems; elem --> 0;) - { - for (mwIndex dim = dims; dim --> 0;) - { - pointData[elem][dim]=data[elem*dims+dim]; - } - } - - mwSize k = 1; - - const mxArray* kData=prhs[2]; - - if (mxIsSparse(kData) || - mxGetNumberOfElements(kData) != 1 || - mxIsComplex(kData) || - !mxIsNumeric(kData)) - { - mexErrMsgIdAndTxt("Damkjer:kannVpTree:prhs", - "First arguement to kannVpTree must be a real valued scalar."); - } - - k = (mwSize)(*((double*)mxGetData(kData))); - - double radius = std::numeric_limits::max(); - - const mxArray* rData=prhs[3]; - - if (mxIsSparse(rData) || - mxGetNumberOfElements(rData) != 1 || - mxIsComplex(rData) || - !mxIsNumeric(rData)) - { - mexErrMsgIdAndTxt("Damkjer:kannVpTree:prhs", - "Second arguement to kannVpTree must be a real valued scalar."); - } - - radius = (*(double*)mxGetData(rData)); - - plhs[0] = mxCreateCellMatrix(elems, 1); - - if (nlhs==2) - { - plhs[1] = mxCreateCellMatrix(elems, 1); - } - -#ifdef _OPENMP - omp_set_dynamic(1); - omp_set_num_threads(omp_get_num_procs()); -#endif - - std::vector, std::deque > > results(pointData.size()); - - #pragma omp parallel for - for (int p = 0; p < pointData.size(); ++p) - { - std::vector q=pointData[p]; - -// std::pair, std::deque > results; - - results[p] = tree.knn(q, k, radius); - } - - #pragma omp parallel for - for (int p = 0; p < pointData.size(); ++p) - { - mwSize neighbors = results[p].first.size(); - - mxArray* neighbor_idxs = 0; - mwIndex* idxs = 0; - mxArray* neighbor_dists = 0; - double* dists = 0; - - #pragma omp critical //(VPSB_KNN_CREATE_IDXS) - { - neighbor_idxs = mxCreateNumericMatrix(0, 0, mxINDEX_CLASS, mxREAL); - mxSetM(neighbor_idxs, neighbors); - mxSetN(neighbor_idxs, 1); - mxSetData(neighbor_idxs, mxMalloc(sizeof(mwIndex)*neighbors*1)); - - idxs = (mwIndex*)mxGetData(neighbor_idxs); - } - - for (mwIndex idx = neighbors; idx --> 0;) - { - idxs[idx]=results[p].first[idx]+1; - } - - #pragma omp critical //(VPSB_KNN_SET_CELL_IDX) - { - mxSetCell(plhs[0], p, neighbor_idxs); - } - - if (nlhs==2) - { - #pragma omp critical //(VPSB_KNN_CREATE_DISTS) - { - neighbor_dists = mxCreateDoubleMatrix(0, 0, mxREAL); - mxSetM(neighbor_dists, neighbors); - mxSetN(neighbor_dists, 1); - mxSetData(neighbor_dists, mxMalloc(sizeof(double)*neighbors*1)); - - dists = mxGetPr(neighbor_dists); - } - - for (mwIndex idx = neighbors; idx --> 0;) - { - dists[idx]=results[p].second[idx]; - } - - #pragma omp critical //(VPSB_KNN_SET_CELL_DISTS) - { - mxSetCell(plhs[1], p, neighbor_dists); - } - } - } -} Index: mkjer/Util/SpatialIndexing/VpTree/newVpTree.cpp =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/newVpTree.cpp (revision 0) +++ (revision ) @@ -1,76 +1,0 @@ -//************************************************************************* -// FILE: newVpTree.cpp -// -// Copyright (C) 2012 Kristian Damkjer. -// -// DESCRIPTION: This MEX source file creates a new VpTree object and -// returns a handle to MATLAB encoded as a raw pointer value. -// -// LIMITATIONS: Unfortunately, the raw pointer trick is necessary to -// persist objects between MEX function calls. Ideally, we -// would simply return the created object. However, that would -// violate MATLAB's current internal memory management system. -// -// SOFTWARE HISTORY: -//> 2012-SEP-11 K. Damkjer -// Initial Coding. -//< -//************************************************************************* - -#ifdef _OPENMP -#include -#endif - -#include "Util/MATLAB/ClassHandle.h" -#include "VpTree.h" - -#ifdef _CHAR16T -#define CHAR16_T -#endif - -#include "mex.h" - -void mexFunction( - int nlhs, mxArray* plhs[], - int nrhs, const mxArray* prhs[]) -{ - if (nrhs < 1 || nrhs > 1) - { - mexErrMsgIdAndTxt("Damkjer:makeVpTree:nargin", - "makeVpTree requires a single input."); - } - - if (nlhs > 1) - { - mexErrMsgIdAndTxt("Damkjer:makeVpTree:nargout", - "makeVpTree requires a single output."); - } - - const mxArray* points=prhs[0]; - - // Check to make sure that points are real numerics - if (mxIsSparse(points) || - mxGetNumberOfDimensions(points) != 2 || - mxIsComplex(points) || - !mxIsNumeric(points)) - { - mexErrMsgIdAndTxt("Damkjer:makeVpTree:prhs", - "Input to makeVpTree must be a full, 2-D matrix representing ND observations."); - } - - const mwSize dims = mxGetM(points); - const mwSize elems = mxGetN(points); - - double* data = mxGetPr(points); - std::deque > pointData(elems, std::vector(dims)); - - for (mwSize elem = elems; elem --> 0;) - { - for (mwSize dim = dims; dim --> 0;) - { - pointData[elem][dim]=data[elem*dims+dim]; - } - } - - plhs[0] = ptrAsMat(new VpTree<>(pointData)); -}