diff --git a/graph/CMakeLists.txt b/felz/CMakeLists.txt similarity index 86% rename from graph/CMakeLists.txt rename to felz/CMakeLists.txt index 260690e..3909905 100644 --- a/graph/CMakeLists.txt +++ b/felz/CMakeLists.txt @@ -23,7 +23,7 @@ include_directories( ${NUMPY_INCLUDE_DIRS} ) -add_library(segment_felz MODULE src/segment.cpp) +add_library(segment_felz MODULE src/segment.cpp src/conversion.cpp) set_target_properties(segment_felz PROPERTIES PREFIX "") target_link_libraries(segment_felz ${Boost_LIBRARIES} @@ -31,7 +31,7 @@ target_link_libraries(segment_felz ${PYTHON_LIBRARIES} ) -add_executable(segment src/segment.cpp) +add_executable(segment src/segment.cpp src/conversion.cpp) target_link_libraries(segment ${Boost_LIBRARIES} ${OpenCV_LIBRARIES} diff --git a/graph/README.md b/felz/README.md similarity index 100% rename from graph/README.md rename to felz/README.md diff --git a/graph/cmake/FindNumPy.cmake b/felz/cmake/FindNumPy.cmake similarity index 100% rename from graph/cmake/FindNumPy.cmake rename to felz/cmake/FindNumPy.cmake diff --git a/graph/src/CVBoostConverter.hpp b/felz/src/CVBoostConverter.hpp similarity index 100% rename from graph/src/CVBoostConverter.hpp rename to felz/src/CVBoostConverter.hpp diff --git a/graph/src/cluster.h b/felz/src/cluster.h similarity index 100% rename from graph/src/cluster.h rename to felz/src/cluster.h diff --git a/felz/src/conversion.cpp b/felz/src/conversion.cpp new file mode 100755 index 0000000..475c670 --- /dev/null +++ b/felz/src/conversion.cpp @@ -0,0 +1,231 @@ +# include "conversion.h" +/* + * The following conversion functions are taken/adapted from OpenCV's cv2.cpp file + * inside modules/python/src2 folder. + */ + +static void init() +{ + import_array(); +} + +static int failmsg(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class PyAllowThreads +{ +public: + PyAllowThreads() : _state(PyEval_SaveThread()) {} + ~PyAllowThreads() + { + PyEval_RestoreThread(_state); + } +private: + PyThreadState* _state; +}; + +class PyEnsureGIL +{ +public: + PyEnsureGIL() : _state(PyGILState_Ensure()) {} + ~PyEnsureGIL() + { + PyGILState_Release(_state); + } +private: + PyGILState_STATE _state; +}; + +using namespace cv; + +static PyObject* failmsgp(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class NumpyAllocator : public MatAllocator +{ +public: + NumpyAllocator() {} + ~NumpyAllocator() {} + + void allocate(int dims, const int* sizes, int type, int*& refcount, + uchar*& datastart, uchar*& data, size_t* step) + { + PyEnsureGIL gil; + + int depth = CV_MAT_DEPTH(type); + int cn = CV_MAT_CN(type); + const int f = (int)(sizeof(size_t)/8); + int typenum = depth == CV_8U ? NPY_UBYTE : depth == CV_8S ? NPY_BYTE : + depth == CV_16U ? NPY_USHORT : depth == CV_16S ? NPY_SHORT : + depth == CV_32S ? NPY_INT : depth == CV_32F ? NPY_FLOAT : + depth == CV_64F ? NPY_DOUBLE : f*NPY_ULONGLONG + (f^1)*NPY_UINT; + int i; + npy_intp _sizes[CV_MAX_DIM+1]; + for( i = 0; i < dims; i++ ) + { + _sizes[i] = sizes[i]; + } + + if( cn > 1 ) + { + _sizes[dims++] = cn; + } + + PyObject* o = PyArray_SimpleNew(dims, _sizes, typenum); + + if(!o) + { + CV_Error_(CV_StsError, ("The numpy array of typenum=%d, ndims=%d can not be created", typenum, dims)); + } + refcount = refcountFromPyObject(o); + + npy_intp* _strides = PyArray_STRIDES(o); + for( i = 0; i < dims - (cn > 1); i++ ) + step[i] = (size_t)_strides[i]; + datastart = data = (uchar*)PyArray_DATA(o); + } + + void deallocate(int* refcount, uchar*, uchar*) + { + PyEnsureGIL gil; + if( !refcount ) + return; + PyObject* o = pyObjectFromRefcount(refcount); + Py_INCREF(o); + Py_DECREF(o); + } +}; + +NumpyAllocator g_numpyAllocator; + +NDArrayConverter::NDArrayConverter() { init(); } + +void NDArrayConverter::init() +{ + import_array(); +} + +cv::Mat NDArrayConverter::toMat(const PyObject *o) +{ + cv::Mat m; + + if(!o || o == Py_None) + { + if( !m.data ) + m.allocator = &g_numpyAllocator; + } + + if( !PyArray_Check(o) ) + { + failmsg("toMat: Object is not a numpy array"); + } + + int typenum = PyArray_TYPE(o); + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT || typenum == NPY_LONG ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if( type < 0 ) + { + failmsg("toMat: Data type = %d is not supported", typenum); + } + + int ndims = PyArray_NDIM(o); + + if(ndims >= CV_MAX_DIM) + { + failmsg("toMat: Dimensionality (=%d) is too high", ndims); + } + + int size[CV_MAX_DIM+1]; + size_t step[CV_MAX_DIM+1], elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(o); + const npy_intp* _strides = PyArray_STRIDES(o); + bool transposed = false; + + for(int i = 0; i < ndims; i++) + { + size[i] = (int)_sizes[i]; + step[i] = (size_t)_strides[i]; + } + + if( ndims == 0 || step[ndims-1] > elemsize ) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if( ndims >= 2 && step[0] < step[1] ) + { + std::swap(size[0], size[1]); + std::swap(step[0], step[1]); + transposed = true; + } + + if( ndims == 3 && size[2] <= CV_CN_MAX && step[1] == elemsize*size[2] ) + { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + + if( ndims > 2) + { + failmsg("toMat: Object has more than 2 dimensions"); + } + + m = Mat(ndims, size, type, PyArray_DATA(o), step); + + if( m.data ) + { + m.refcount = refcountFromPyObject(o); + m.addref(); // protect the original numpy array from deallocation + // (since Mat destructor will decrement the reference counter) + }; + m.allocator = &g_numpyAllocator; + + if( transposed ) + { + Mat tmp; + tmp.allocator = &g_numpyAllocator; + transpose(m, tmp); + m = tmp; + } + return m; +} + +PyObject* NDArrayConverter::toNDArray(const cv::Mat& m) +{ + if( !m.data ) + Py_RETURN_NONE; + Mat temp, *p = (Mat*)&m; + if(!p->refcount || p->allocator != &g_numpyAllocator) + { + temp.allocator = &g_numpyAllocator; + m.copyTo(temp); + p = &temp; + } + p->addref(); + return pyObjectFromRefcount(p->refcount); +} diff --git a/felz/src/conversion.h b/felz/src/conversion.h new file mode 100755 index 0000000..3975e17 --- /dev/null +++ b/felz/src/conversion.h @@ -0,0 +1,60 @@ +# ifndef __COVERSION_OPENCV_H__ +# define __COVERSION_OPENCV_H__ + +#include +#include +#include +#include +#include "numpy/ndarrayobject.h" + +static PyObject* opencv_error = 0; + +static int failmsg(const char *fmt, ...); + +class PyAllowThreads; + +class PyEnsureGIL; + +#define ERRWRAP2(expr) \ +try \ +{ \ + PyAllowThreads allowThreads; \ + expr; \ +} \ +catch (const cv::Exception &e) \ +{ \ + PyErr_SetString(opencv_error, e.what()); \ + return 0; \ +} + +static PyObject* failmsgp(const char *fmt, ...); + +static size_t REFCOUNT_OFFSET = (size_t)&(((PyObject*)0)->ob_refcnt) + + (0x12345678 != *(const size_t*)"\x78\x56\x34\x12\0\0\0\0\0")*sizeof(int); + +static inline PyObject* pyObjectFromRefcount(const int* refcount) +{ + return (PyObject*)((size_t)refcount - REFCOUNT_OFFSET); +} + +static inline int* refcountFromPyObject(const PyObject* obj) +{ + return (int*)((size_t)obj + REFCOUNT_OFFSET); +} + + +class NumpyAllocator; + +enum { ARG_NONE = 0, ARG_MAT = 1, ARG_SCALAR = 2 }; + +class NDArrayConverter +{ +private: + void init(); +public: + NDArrayConverter(); + cv::Mat toMat(const PyObject* o); + PyObject* toNDArray(const cv::Mat& mat); +}; + +# endif diff --git a/graph/src/edge.h b/felz/src/edge.h similarity index 100% rename from graph/src/edge.h rename to felz/src/edge.h diff --git a/graph/src/segment-graph.h b/felz/src/segment-graph.h similarity index 100% rename from graph/src/segment-graph.h rename to felz/src/segment-graph.h diff --git a/graph/src/segment-image.h b/felz/src/segment-image.h similarity index 86% rename from graph/src/segment-image.h rename to felz/src/segment-image.h index 588cd26..5d03233 100644 --- a/graph/src/segment-image.h +++ b/felz/src/segment-image.h @@ -5,7 +5,7 @@ #include #define square(x) ((x)*(x)) -#define COLORED_OUTPUT +//#define COLORED_OUTPUT // random color cv::Vec3b random_rgb(){ @@ -19,10 +19,12 @@ cv::Vec3b random_rgb(){ /// Difference between two pixels p1, p2 in image static inline float diff(cv::Mat image, int p1, int p2) { - //cv::Vec3b d = image.at(p1) - image.at(p2); - cv::Vec3f a = image.at(p1); - cv::Vec3f b = image.at(p2); - return sqrt(square(a[0] - b[0]) + square(a[1] - b[1]) + square(a[2] - b[2])); + if (image.channels() == 3) { + cv::Vec3f a = image.at(p1); + cv::Vec3f b = image.at(p2); + return sqrt(square(a[0] - b[0]) + square(a[1] - b[1]) + square(a[2] - b[2])); + } else + return fabs(image.at(p1) - image.at(p2)); } /* @@ -36,7 +38,7 @@ static inline float diff(cv::Mat image, int p1, int p2) { * min_size: minimum component size (enforced by post-processing stage). */ cv::Mat segment_image(cv::Mat image, float sigma, float c, int min_size) { - image.convertTo(image, CV_32FC3); + image.convertTo(image, CV_MAKETYPE(CV_32F, image.channels())); cv::GaussianBlur(image, image, cv::Size(5, 5), sigma); //image = smooth(image, sigma); diff --git a/graph/src/segment.cpp b/felz/src/segment.cpp similarity index 57% rename from graph/src/segment.cpp rename to felz/src/segment.cpp index 718166e..6cc585c 100644 --- a/graph/src/segment.cpp +++ b/felz/src/segment.cpp @@ -4,13 +4,16 @@ #include #include "segment-image.h" #include -#include "CVBoostConverter.hpp" +//#include "CVBoostConverter.hpp" +#include "conversion.h" #include using namespace boost::python; -cv::Mat segment(cv::Mat image, float sigma, int k, int min_size) { - return segment_image(image, sigma, k, min_size); +PyObject * segment(PyObject * image_, float sigma, int k, int min_size) { + NDArrayConverter cvt; + cv::Mat image = cvt.toMat(image_); + return cvt.toNDArray(segment_image(image, sigma, k, min_size)); } static void init_ar() { @@ -22,16 +25,16 @@ BOOST_PYTHON_MODULE(segment_felz) { init_ar(); //initialize converters - to_python_converter(); - bcvt::matFromNDArrayBoostConverter(); + bcvt::matFromNDArrayBoostConverter();*/ def("segment", segment); } int main(int argc, char **argv) { if (argc != 6) { - std::cerr << "usage: " << argv[0] << " sigma k min input output" << std::endl; + std::cerr << "usage: " << argv[0] << " sigma k min image output" << std::endl; return 1; } @@ -39,18 +42,18 @@ int main(int argc, char **argv) { int k = atoi(argv[2]); int min_size = atoi(argv[3]); - cv::Mat input = cv::imread(argv[4]); - //cv::resize(input, input, cv::Size(input.cols * 0.4, input.rows * 0.4)); + cv::Mat image = cv::imread(argv[4]); + //cv::resize(image, image, cv::Size(image.cols * 0.4, image.rows * 0.4)); //int num_ccs; std::cout << "Processing..." << std::endl; std::clock_t begin = std::clock(); - cv::Mat seg = segment(input, sigma, k, min_size); + cv::Mat seg = segment_image(image, sigma, k, min_size); std::clock_t end = std::clock(); double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC; cv::imwrite(argv[5], seg); - std::cout << "Times passed in seconds: " << elapsed_secs << std::endl; + std::cout << "Time passed in seconds: " << elapsed_secs << std::endl; //std::cout << "Got " << num_ccs << " components." << std::endl; return 0; diff --git a/random/CMakeLists.txt b/random/CMakeLists.txt new file mode 100644 index 0000000..dda7375 --- /dev/null +++ b/random/CMakeLists.txt @@ -0,0 +1,40 @@ +cmake_minimum_required(VERSION 2.6) +project(randomselection) +set(PROJECT_NAME randomselection) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +# Boost +set(Boost_USE_STATIC_LIBS ON) +set(Boost_USE_MULTITHREADED ON) +find_package(Boost COMPONENTS python REQUIRED) + +find_package(OpenCV REQUIRED) +find_package(PythonLibs REQUIRED) +find_package(NumPy REQUIRED) +find_package(Objectness REQUIRED) + +set(CMAKE_CXX_FLAGS "-std=c++11 -g") + +include_directories( + include + ${PYTHON_INCLUDE_DIRS} + ${NUMPY_INCLUDE_DIRS} +) + +add_library(random_selection MODULE src/main.cpp src/conversion.cpp) +set_target_properties(random_selection PROPERTIES PREFIX "") +target_link_libraries(random_selection + ${Boost_LIBRARIES} + ${OpenCV_LIBRARIES} + ${PYTHON_LIBRARIES} + ${OBJECTNESS_LIBRARIES} +) + +add_executable(random src/main.cpp src/conversion.cpp) +target_link_libraries(random + ${Boost_LIBRARIES} + ${OpenCV_LIBRARIES} + ${PYTHON_LIBRARIES} + ${OBJECTNESS_LIBRARIES} +) diff --git a/random/cmake/FindNumPy.cmake b/random/cmake/FindNumPy.cmake new file mode 100644 index 0000000..eafed16 --- /dev/null +++ b/random/cmake/FindNumPy.cmake @@ -0,0 +1,102 @@ +# - Find the NumPy libraries +# This module finds if NumPy is installed, and sets the following variables +# indicating where it is. +# +# TODO: Update to provide the libraries and paths for linking npymath lib. +# +# NUMPY_FOUND - was NumPy found +# NUMPY_VERSION - the version of NumPy found as a string +# NUMPY_VERSION_MAJOR - the major version number of NumPy +# NUMPY_VERSION_MINOR - the minor version number of NumPy +# NUMPY_VERSION_PATCH - the patch version number of NumPy +# NUMPY_VERSION_DECIMAL - e.g. version 1.6.1 is 10601 +# NUMPY_INCLUDE_DIRS - path to the NumPy include files + +#============================================================================ +# Copyright 2012 Continuum Analytics, Inc. +# +# MIT License +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files +# (the "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to permit +# persons to whom the Software is furnished to do so, subject to +# the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR +# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +# OTHER DEALINGS IN THE SOFTWARE. +# +#============================================================================ + +# Finding NumPy involves calling the Python interpreter +if(NumPy_FIND_REQUIRED) + find_package(PythonInterp REQUIRED) +else() + find_package(PythonInterp) +endif() + +if(NOT PYTHONINTERP_FOUND) + set(NUMPY_FOUND FALSE) + return() +endif() + +execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-c" + "import numpy as n; print(n.__version__); print(n.get_include());" + RESULT_VARIABLE _NUMPY_SEARCH_SUCCESS + OUTPUT_VARIABLE _NUMPY_VALUES_OUTPUT + ERROR_VARIABLE _NUMPY_ERROR_VALUE + OUTPUT_STRIP_TRAILING_WHITESPACE) + +if(NOT _NUMPY_SEARCH_SUCCESS MATCHES 0) + if(NumPy_FIND_REQUIRED) + message(FATAL_ERROR + "NumPy import failure:\n${_NUMPY_ERROR_VALUE}") + endif() + set(NUMPY_FOUND FALSE) + return() +endif() + +# Convert the process output into a list +string(REGEX REPLACE ";" "\\\\;" _NUMPY_VALUES ${_NUMPY_VALUES_OUTPUT}) +string(REGEX REPLACE "\n" ";" _NUMPY_VALUES ${_NUMPY_VALUES}) +# Just in case there is unexpected output from the Python command. +list(GET _NUMPY_VALUES -2 NUMPY_VERSION) +list(GET _NUMPY_VALUES -1 NUMPY_INCLUDE_DIRS) + +string(REGEX MATCH "^[0-9]+\\.[0-9]+\\.[0-9]+" _VER_CHECK "${NUMPY_VERSION}") +if("${_VER_CHECK}" STREQUAL "") + # The output from Python was unexpected. Raise an error always + # here, because we found NumPy, but it appears to be corrupted somehow. + message(FATAL_ERROR + "Requested version and include path from NumPy, got instead:\n${_NUMPY_VALUES_OUTPUT}\n") + return() +endif() + +# Make sure all directory separators are '/' +string(REGEX REPLACE "\\\\" "/" NUMPY_INCLUDE_DIRS ${NUMPY_INCLUDE_DIRS}) + +# Get the major and minor version numbers +string(REGEX REPLACE "\\." ";" _NUMPY_VERSION_LIST ${NUMPY_VERSION}) +list(GET _NUMPY_VERSION_LIST 0 NUMPY_VERSION_MAJOR) +list(GET _NUMPY_VERSION_LIST 1 NUMPY_VERSION_MINOR) +list(GET _NUMPY_VERSION_LIST 2 NUMPY_VERSION_PATCH) +string(REGEX MATCH "[0-9]*" NUMPY_VERSION_PATCH ${NUMPY_VERSION_PATCH}) +math(EXPR NUMPY_VERSION_DECIMAL + "(${NUMPY_VERSION_MAJOR} * 10000) + (${NUMPY_VERSION_MINOR} * 100) + ${NUMPY_VERSION_PATCH}") + +find_package_message(NUMPY + "Found NumPy: version \"${NUMPY_VERSION}\" ${NUMPY_INCLUDE_DIRS}" + "${NUMPY_INCLUDE_DIRS}${NUMPY_VERSION}") + +set(NUMPY_FOUND TRUE) + diff --git a/random/cmake/FindObjectness.cmake b/random/cmake/FindObjectness.cmake new file mode 100644 index 0000000..a6e5241 --- /dev/null +++ b/random/cmake/FindObjectness.cmake @@ -0,0 +1,24 @@ +# - Try to find Objectness +# Once done this will define +# OBJECTNESS_FOUND - System has Objectness +# OBJECTNESS_INCLUDE_DIRS - The Objectness include directories +# OBJECTNESS_LIBRARIES - The libraries needed to use LibXml2 + +find_path(OBJECTNESS_INCLUDE_DIR Objectness/Objectness.h + PATHS ${CMAKE_CURRENT_SOURCE_DIR}/include +) + +find_library(OBJECTNESS_LIBRARY NAMES objectness + PATHS ${CMAKE_CURRENT_SOURCE_DIR}/lib +) + +set(OBJECTNESS_LIBRARIES ${OBJECTNESS_LIBRARY} ) +set(OBJECTNESS_INCLUDE_DIRS ${OBJECTNESS_INCLUDE_DIR} ) + +if (OBJECTNESS_INCLUDE_DIR AND OBJECTNESS_LIBRARIES) + set(OBJECTNESS_FOUND true) +else() + set(OBJECTNESS_FOUND false) +endif() + +mark_as_advanced(OBJECTNESS_FOUND) diff --git a/random/include/LibLinear/CMakeLists.txt b/random/include/LibLinear/CMakeLists.txt new file mode 100644 index 0000000..a63d14e --- /dev/null +++ b/random/include/LibLinear/CMakeLists.txt @@ -0,0 +1,6 @@ +#Do not include files (from libraries) with a main function to avoid conflict +#with the main function in the project. +set(LIBLINEAR_LIST linear.cpp tron.cpp) +set(BLAS_LIST blas/daxpy.c blas/ddot.c blas/dnrm2.c blas/dscal.c) +add_library(LIBLINEAR ${LIBLINEAR_LIST}) +add_library(BLAS ${BLAS_LIST}) diff --git a/random/include/LibLinear/LibLinear.vcxproj b/random/include/LibLinear/LibLinear.vcxproj new file mode 100755 index 0000000..62d1e65 --- /dev/null +++ b/random/include/LibLinear/LibLinear.vcxproj @@ -0,0 +1,94 @@ + + + + + Debug + x64 + + + Release + x64 + + + + {86266F16-8B7E-4666-A12F-96E351579ADA} + Win32Proj + LibLinear + + + + StaticLibrary + true + MultiByte + v110 + + + StaticLibrary + false + true + MultiByte + v110 + + + + + + + + + + + + + + + + + Level3 + Disabled + WIN32;_DEBUG;_LIB;%(PreprocessorDefinitions) + + + Console + true + + + + + Level3 + + + MaxSpeed + true + true + WIN32;NDEBUG;_LIB;%(PreprocessorDefinitions) + + + Console + true + true + true + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/random/include/LibLinear/LibLinear.vcxproj.filters b/random/include/LibLinear/LibLinear.vcxproj.filters new file mode 100755 index 0000000..2c20724 --- /dev/null +++ b/random/include/LibLinear/LibLinear.vcxproj.filters @@ -0,0 +1,57 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms + + + + + Source Files + + + Source Files + + + Source Files + + + BALAS + + + BALAS + + + BALAS + + + BALAS + + + + + Header Files + + + Header Files + + + BALAS + + + BALAS + + + + + + \ No newline at end of file diff --git a/random/include/LibLinear/README.1.93.txt b/random/include/LibLinear/README.1.93.txt new file mode 100755 index 0000000..9b5c039 --- /dev/null +++ b/random/include/LibLinear/README.1.93.txt @@ -0,0 +1,531 @@ +LIBLINEAR is a simple package for solving large-scale regularized linear +classification and regression. It currently supports +- L2-regularized logistic regression/L2-loss support vector classification/L1-loss support vector classification +- L1-regularized L2-loss support vector classification/L1-regularized logistic regression +- L2-regularized L2-loss support vector regression/L1-loss support vector regression. +This document explains the usage of LIBLINEAR. + +To get started, please read the ``Quick Start'' section first. +For developers, please check the ``Library Usage'' section to learn +how to integrate LIBLINEAR in your software. + +Table of Contents +================= + +- When to use LIBLINEAR but not LIBSVM +- Quick Start +- Installation +- `train' Usage +- `predict' Usage +- Examples +- Library Usage +- Building Windows Binaries +- Additional Information +- MATLAB/OCTAVE interface +- PYTHON interface + +When to use LIBLINEAR but not LIBSVM +==================================== + +There are some large data for which with/without nonlinear mappings +gives similar performances. Without using kernels, one can +efficiently train a much larger set via linear classification/regression. +These data usually have a large number of features. Document classification +is an example. + +Warning: While generally liblinear is very fast, its default solver +may be slow under certain situations (e.g., data not scaled or C is +large). See Appendix B of our SVM guide about how to handle such +cases. +http://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf + +Warning: If you are a beginner and your data sets are not large, you +should consider LIBSVM first. + +LIBSVM page: +http://www.csie.ntu.edu.tw/~cjlin/libsvm + + +Quick Start +=========== + +See the section ``Installation'' for installing LIBLINEAR. + +After installation, there are programs `train' and `predict' for +training and testing, respectively. + +About the data format, please check the README file of LIBSVM. Note +that feature index must start from 1 (but not 0). + +A sample classification data included in this package is `heart_scale'. + +Type `train heart_scale', and the program will read the training +data and output the model file `heart_scale.model'. If you have a test +set called heart_scale.t, then type `predict heart_scale.t +heart_scale.model output' to see the prediction accuracy. The `output' +file contains the predicted class labels. + +For more information about `train' and `predict', see the sections +`train' Usage and `predict' Usage. + +To obtain good performances, sometimes one needs to scale the +data. Please check the program `svm-scale' of LIBSVM. For large and +sparse data, use `-l 0' to keep the sparsity. + +Installation +============ + +On Unix systems, type `make' to build the `train' and `predict' +programs. Run them without arguments to show the usages. + +On other systems, consult `Makefile' to build them (e.g., see +'Building Windows binaries' in this file) or use the pre-built +binaries (Windows binaries are in the directory `windows'). + +This software uses some level-1 BLAS subroutines. The needed functions are +included in this package. If a BLAS library is available on your +machine, you may use it by modifying the Makefile: Unmark the following line + + #LIBS ?= -lblas + +and mark + + LIBS ?= blas/blas.a + +`train' Usage +============= + +Usage: train [options] training_set_file [model_file] +options: +-s type : set type of solver (default 1) + for multi-class classification + 0 -- L2-regularized logistic regression (primal) + 1 -- L2-regularized L2-loss support vector classification (dual) + 2 -- L2-regularized L2-loss support vector classification (primal) + 3 -- L2-regularized L1-loss support vector classification (dual) + 4 -- support vector classification by Crammer and Singer + 5 -- L1-regularized L2-loss support vector classification + 6 -- L1-regularized logistic regression + 7 -- L2-regularized logistic regression (dual) + for regression + 11 -- L2-regularized L2-loss support vector regression (primal) + 12 -- L2-regularized L2-loss support vector regression (dual) + 13 -- L2-regularized L1-loss support vector regression (dual) +-c cost : set the parameter C (default 1) +-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1) +-e epsilon : set tolerance of termination criterion + -s 0 and 2 + |f'(w)|_2 <= eps*min(pos,neg)/l*|f'(w0)|_2, + where f is the primal function and pos/neg are # of + positive/negative data (default 0.01) + -s 11 + |f'(w)|_2 <= eps*|f'(w0)|_2 (default 0.001) + -s 1, 3, 4 and 7 + Dual maximal violation <= eps; similar to libsvm (default 0.1) + -s 5 and 6 + |f'(w)|_inf <= eps*min(pos,neg)/l*|f'(w0)|_inf, + where f is the primal function (default 0.01) + -s 12 and 13\n" + |f'(alpha)|_1 <= eps |f'(alpha0)|, + where f is the dual function (default 0.1) +-B bias : if bias >= 0, instance x becomes [x; bias]; if < 0, no bias term added (default -1) +-wi weight: weights adjust the parameter C of different classes (see README for details) +-v n: n-fold cross validation mode +-q : quiet mode (no outputs) + +Option -v randomly splits the data into n parts and calculates cross +validation accuracy on them. + +Formulations: + +For L2-regularized logistic regression (-s 0), we solve + +min_w w^Tw/2 + C \sum log(1 + exp(-y_i w^Tx_i)) + +For L2-regularized L2-loss SVC dual (-s 1), we solve + +min_alpha 0.5(alpha^T (Q + I/2/C) alpha) - e^T alpha + s.t. 0 <= alpha_i, + +For L2-regularized L2-loss SVC (-s 2), we solve + +min_w w^Tw/2 + C \sum max(0, 1- y_i w^Tx_i)^2 + +For L2-regularized L1-loss SVC dual (-s 3), we solve + +min_alpha 0.5(alpha^T Q alpha) - e^T alpha + s.t. 0 <= alpha_i <= C, + +For L1-regularized L2-loss SVC (-s 5), we solve + +min_w \sum |w_j| + C \sum max(0, 1- y_i w^Tx_i)^2 + +For L1-regularized logistic regression (-s 6), we solve + +min_w \sum |w_j| + C \sum log(1 + exp(-y_i w^Tx_i)) + +For L2-regularized logistic regression (-s 7), we solve + +min_alpha 0.5(alpha^T Q alpha) + \sum alpha_i*log(alpha_i) + \sum (C-alpha_i)*log(C-alpha_i) - a constant + s.t. 0 <= alpha_i <= C, + +where + +Q is a matrix with Q_ij = y_i y_j x_i^T x_j. + +For L2-regularized L2-loss SVR (-s 11), we solve + +min_w w^Tw/2 + C \sum max(0, |y_i-w^Tx_i|-epsilon)^2 + +For L2-regularized L2-loss SVR dual (-s 12), we solve + +min_beta 0.5(beta^T (Q + lambda I/2/C) beta) - y^T beta + \sum |beta_i| + +For L2-regularized L1-loss SVR dual (-s 13), we solve + +min_beta 0.5(beta^T Q beta) - y^T beta + \sum |beta_i| + s.t. -C <= beta_i <= C, + +where + +Q is a matrix with Q_ij = x_i^T x_j. + +If bias >= 0, w becomes [w; w_{n+1}] and x becomes [x; bias]. + +The primal-dual relationship implies that -s 1 and -s 2 give the same +model, -s 0 and -s 7 give the same, and -s 11 and -s 12 give the same. + +We implement 1-vs-the rest multi-class strategy for classification. +In training i vs. non_i, their C parameters are (weight from -wi)*C +and C, respectively. If there are only two classes, we train only one +model. Thus weight1*C vs. weight2*C is used. See examples below. + +We also implement multi-class SVM by Crammer and Singer (-s 4): + +min_{w_m, \xi_i} 0.5 \sum_m ||w_m||^2 + C \sum_i \xi_i + s.t. w^T_{y_i} x_i - w^T_m x_i >= \e^m_i - \xi_i \forall m,i + +where e^m_i = 0 if y_i = m, + e^m_i = 1 if y_i != m, + +Here we solve the dual problem: + +min_{\alpha} 0.5 \sum_m ||w_m(\alpha)||^2 + \sum_i \sum_m e^m_i alpha^m_i + s.t. \alpha^m_i <= C^m_i \forall m,i , \sum_m \alpha^m_i=0 \forall i + +where w_m(\alpha) = \sum_i \alpha^m_i x_i, +and C^m_i = C if m = y_i, + C^m_i = 0 if m != y_i. + +`predict' Usage +=============== + +Usage: predict [options] test_file model_file output_file +options: +-b probability_estimates: whether to output probability estimates, 0 or 1 (default 0); currently for logistic regression only +-q : quiet mode (no outputs) + +Note that -b is only needed in the prediction phase. This is different +from the setting of LIBSVM. + +Examples +======== + +> train data_file + +Train linear SVM with L2-loss function. + +> train -s 0 data_file + +Train a logistic regression model. + +> train -v 5 -e 0.001 data_file + +Do five-fold cross-validation using L2-loss svm. +Use a smaller stopping tolerance 0.001 than the default +0.1 if you want more accurate solutions. + +> train -c 10 -w1 2 -w2 5 -w3 2 four_class_data_file + +Train four classifiers: +positive negative Cp Cn +class 1 class 2,3,4. 20 10 +class 2 class 1,3,4. 50 10 +class 3 class 1,2,4. 20 10 +class 4 class 1,2,3. 10 10 + +> train -c 10 -w3 1 -w2 5 two_class_data_file + +If there are only two classes, we train ONE model. +The C values for the two classes are 10 and 50. + +> predict -b 1 test_file data_file.model output_file + +Output probability estimates (for logistic regression only). + +Library Usage +============= + +- Function: model* train(const struct problem *prob, + const struct parameter *param); + + This function constructs and returns a linear classification + or regression model according to the given training data and + parameters. + + struct problem describes the problem: + + struct problem + { + int l, n; + int *y; + struct feature_node **x; + double bias; + }; + + where `l' is the number of training data. If bias >= 0, we assume + that one additional feature is added to the end of each data + instance. `n' is the number of feature (including the bias feature + if bias >= 0). `y' is an array containing the target values. (integers + in classification, real numbers in regression) And `x' is an array + of pointers, each of which points to a sparse representation (array + of feature_node) of one training vector. + + For example, if we have the following training data: + + LABEL ATTR1 ATTR2 ATTR3 ATTR4 ATTR5 + ----- ----- ----- ----- ----- ----- + 1 0 0.1 0.2 0 0 + 2 0 0.1 0.3 -1.2 0 + 1 0.4 0 0 0 0 + 2 0 0.1 0 1.4 0.5 + 3 -0.1 -0.2 0.1 1.1 0.1 + + and bias = 1, then the components of problem are: + + l = 5 + n = 6 + + y -> 1 2 1 2 3 + + x -> [ ] -> (2,0.1) (3,0.2) (6,1) (-1,?) + [ ] -> (2,0.1) (3,0.3) (4,-1.2) (6,1) (-1,?) + [ ] -> (1,0.4) (6,1) (-1,?) + [ ] -> (2,0.1) (4,1.4) (5,0.5) (6,1) (-1,?) + [ ] -> (1,-0.1) (2,-0.2) (3,0.1) (4,1.1) (5,0.1) (6,1) (-1,?) + + struct parameter describes the parameters of a linear classification + or regression model: + + struct parameter + { + int solver_type; + + /* these are for training only */ + double eps; /* stopping criteria */ + double C; + int nr_weight; + int *weight_label; + double* weight; + double p; + }; + + solver_type can be one of L2R_LR, L2R_L2LOSS_SVC_DUAL, L2R_L2LOSS_SVC, L2R_L1LOSS_SVC_DUAL, MCSVM_CS, L1R_L2LOSS_SVC, L1R_LR, L2R_LR_DUAL, L2R_L2LOSS_SVR, L2R_L2LOSS_SVR_DUAL, L2R_L1LOSS_SVR_DUAL. + for classification + L2R_LR L2-regularized logistic regression (primal) + L2R_L2LOSS_SVC_DUAL L2-regularized L2-loss support vector classification (dual) + L2R_L2LOSS_SVC L2-regularized L2-loss support vector classification (primal) + L2R_L1LOSS_SVC_DUAL L2-regularized L1-loss support vector classification (dual) + MCSVM_CS support vector classification by Crammer and Singer + L1R_L2LOSS_SVC L1-regularized L2-loss support vector classification + L1R_LR L1-regularized logistic regression + L2R_LR_DUAL L2-regularized logistic regression (dual) + for regression + L2R_L2LOSS_SVR L2-regularized L2-loss support vector regression (primal) + L2R_L2LOSS_SVR_DUAL L2-regularized L2-loss support vector regression (dual) + L2R_L1LOSS_SVR_DUAL L2-regularized L1-loss support vector regression (dual) + + C is the cost of constraints violation. + p is the sensitiveness of loss of support vector regression. + eps is the stopping criterion. + + nr_weight, weight_label, and weight are used to change the penalty + for some classes (If the weight for a class is not changed, it is + set to 1). This is useful for training classifier using unbalanced + input data or with asymmetric misclassification cost. + + nr_weight is the number of elements in the array weight_label and + weight. Each weight[i] corresponds to weight_label[i], meaning that + the penalty of class weight_label[i] is scaled by a factor of weight[i]. + + If you do not want to change penalty for any of the classes, + just set nr_weight to 0. + + *NOTE* To avoid wrong parameters, check_parameter() should be + called before train(). + + struct model stores the model obtained from the training procedure: + + struct model + { + struct parameter param; + int nr_class; /* number of classes */ + int nr_feature; + double *w; + int *label; /* label of each class */ + double bias; + }; + + param describes the parameters used to obtain the model. + + nr_class and nr_feature are the number of classes and features, + respectively. nr_class = 2 for regression. + + The nr_feature*nr_class array w gives feature weights. We use one + against the rest for multi-class classification, so each feature + index corresponds to nr_class weight values. Weights are + organized in the following way + + +------------------+------------------+------------+ + | nr_class weights | nr_class weights | ... + | for 1st feature | for 2nd feature | + +------------------+------------------+------------+ + + If bias >= 0, x becomes [x; bias]. The number of features is + increased by one, so w is a (nr_feature+1)*nr_class array. The + value of bias is stored in the variable bias. + + The array label stores class labels. + +- Function: void cross_validation(const problem *prob, const parameter *param, int nr_fold, double *target); + + This function conducts cross validation. Data are separated to + nr_fold folds. Under given parameters, sequentially each fold is + validated using the model from training the remaining. Predicted + labels in the validation process are stored in the array called + target. + + The format of prob is same as that for train(). + +- Function: double predict(const model *model_, const feature_node *x); + + For a classification model, the predicted class for x is returned. + For a regression model, the function value of x calculated using + the model is returned. + +- Function: double predict_values(const struct model *model_, + const struct feature_node *x, double* dec_values); + + This function gives nr_w decision values in the array dec_values. + nr_w=1 if regression is applied or the number of classes is two. An exception is + multi-class svm by Crammer and Singer (-s 4), where nr_w = 2 if there are two classes. For all other situations, nr_w is the + number of classes. + + We implement one-vs-the rest multi-class strategy (-s 0,1,2,3,5,6,7) + and multi-class svm by Crammer and Singer (-s 4) for multi-class SVM. + The class with the highest decision value is returned. + +- Function: double predict_probability(const struct model *model_, + const struct feature_node *x, double* prob_estimates); + + This function gives nr_class probability estimates in the array + prob_estimates. nr_class can be obtained from the function + get_nr_class. The class with the highest probability is + returned. Currently, we support only the probability outputs of + logistic regression. + +- Function: int get_nr_feature(const model *model_); + + The function gives the number of attributes of the model. + +- Function: int get_nr_class(const model *model_); + + The function gives the number of classes of the model. + For a regression model, 2 is returned. + +- Function: void get_labels(const model *model_, int* label); + + This function outputs the name of labels into an array called label. + For a regression model, label is unchanged. + +- Function: const char *check_parameter(const struct problem *prob, + const struct parameter *param); + + This function checks whether the parameters are within the feasible + range of the problem. This function should be called before calling + train() and cross_validation(). It returns NULL if the + parameters are feasible, otherwise an error message is returned. + +- Function: int save_model(const char *model_file_name, + const struct model *model_); + + This function saves a model to a file; returns 0 on success, or -1 + if an error occurs. + +- Function: struct model *load_model(const char *model_file_name); + + This function returns a pointer to the model read from the file, + or a null pointer if the model could not be loaded. + +- Function: void free_model_content(struct model *model_ptr); + + This function frees the memory used by the entries in a model structure. + +- Function: void free_and_destroy_model(struct model **model_ptr_ptr); + + This function frees the memory used by a model and destroys the model + structure. + +- Function: void destroy_param(struct parameter *param); + + This function frees the memory used by a parameter set. + +- Function: void set_print_string_function(void (*print_func)(const char *)); + + Users can specify their output format by a function. Use + set_print_string_function(NULL); + for default printing to stdout. + +Building Windows Binaries +========================= + +Windows binaries are in the directory `windows'. To build them via +Visual C++, use the following steps: + +1. Open a dos command box and change to liblinear directory. If +environment variables of VC++ have not been set, type + +"C:\Program Files\Microsoft Visual Studio 10.0\VC\bin\vcvars32.bat" + +You may have to modify the above command according which version of +VC++ or where it is installed. + +2. Type + +nmake -f Makefile.win clean all + + +MATLAB/OCTAVE Interface +======================= + +Please check the file README in the directory `matlab'. + +PYTHON Interface +================ + +Please check the file README in the directory `python'. + +Additional Information +====================== + +If you find LIBLINEAR helpful, please cite it as + +R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang, and C.-J. Lin. +LIBLINEAR: A Library for Large Linear Classification, Journal of +Machine Learning Research 9(2008), 1871-1874. Software available at +http://www.csie.ntu.edu.tw/~cjlin/liblinear + +For any questions and comments, please send your email to +cjlin@csie.ntu.edu.tw + + diff --git a/random/include/LibLinear/blas/blas.h b/random/include/LibLinear/blas/blas.h new file mode 100755 index 0000000..67b7009 --- /dev/null +++ b/random/include/LibLinear/blas/blas.h @@ -0,0 +1,25 @@ +/* blas.h -- C header file for BLAS Ver 1.0 */ +/* Jesse Bennett March 23, 2000 */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef BLAS_INCLUDE +#define BLAS_INCLUDE + +/* Data types specific to BLAS implementation */ +typedef struct { float r, i; } fcomplex; +typedef struct { double r, i; } dcomplex; +typedef int blasbool; + +#include "blasp.h" /* Prototypes for all BLAS functions */ + +#define FALSE 0 +#define TRUE 1 + +/* Macro functions */ +#define MIN(a,b) ((a) <= (b) ? (a) : (b)) +#define MAX(a,b) ((a) >= (b) ? (a) : (b)) + +#endif diff --git a/random/include/LibLinear/blas/blasp.h b/random/include/LibLinear/blas/blasp.h new file mode 100755 index 0000000..ff244d8 --- /dev/null +++ b/random/include/LibLinear/blas/blasp.h @@ -0,0 +1,430 @@ +/* blasp.h -- C prototypes for BLAS Ver 1.0 */ +/* Jesse Bennett March 23, 2000 */ + +/* Functions listed in alphabetical order */ + +#ifdef F2C_COMPAT + +void cdotc_(fcomplex *dotval, int *n, fcomplex *cx, int *incx, + fcomplex *cy, int *incy); + +void cdotu_(fcomplex *dotval, int *n, fcomplex *cx, int *incx, + fcomplex *cy, int *incy); + +double sasum_(int *n, float *sx, int *incx); + +double scasum_(int *n, fcomplex *cx, int *incx); + +double scnrm2_(int *n, fcomplex *x, int *incx); + +double sdot_(int *n, float *sx, int *incx, float *sy, int *incy); + +double snrm2_(int *n, float *x, int *incx); + +void zdotc_(dcomplex *dotval, int *n, dcomplex *cx, int *incx, + dcomplex *cy, int *incy); + +void zdotu_(dcomplex *dotval, int *n, dcomplex *cx, int *incx, + dcomplex *cy, int *incy); + +#else + +fcomplex cdotc_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); + +fcomplex cdotu_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); + +float sasum_(int *n, float *sx, int *incx); + +float scasum_(int *n, fcomplex *cx, int *incx); + +float scnrm2_(int *n, fcomplex *x, int *incx); + +float sdot_(int *n, float *sx, int *incx, float *sy, int *incy); + +float snrm2_(int *n, float *x, int *incx); + +dcomplex zdotc_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); + +dcomplex zdotu_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); + +#endif + +/* Remaining functions listed in alphabetical order */ + +int caxpy_(int *n, fcomplex *ca, fcomplex *cx, int *incx, fcomplex *cy, + int *incy); + +int ccopy_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); + +int cgbmv_(char *trans, int *m, int *n, int *kl, int *ku, + fcomplex *alpha, fcomplex *a, int *lda, fcomplex *x, int *incx, + fcomplex *beta, fcomplex *y, int *incy); + +int cgemm_(char *transa, char *transb, int *m, int *n, int *k, + fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, int *ldb, + fcomplex *beta, fcomplex *c, int *ldc); + +int cgemv_(char *trans, int *m, int *n, fcomplex *alpha, fcomplex *a, + int *lda, fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, + int *incy); + +int cgerc_(int *m, int *n, fcomplex *alpha, fcomplex *x, int *incx, + fcomplex *y, int *incy, fcomplex *a, int *lda); + +int cgeru_(int *m, int *n, fcomplex *alpha, fcomplex *x, int *incx, + fcomplex *y, int *incy, fcomplex *a, int *lda); + +int chbmv_(char *uplo, int *n, int *k, fcomplex *alpha, fcomplex *a, + int *lda, fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, + int *incy); + +int chemm_(char *side, char *uplo, int *m, int *n, fcomplex *alpha, + fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, + fcomplex *c, int *ldc); + +int chemv_(char *uplo, int *n, fcomplex *alpha, fcomplex *a, int *lda, + fcomplex *x, int *incx, fcomplex *beta, fcomplex *y, int *incy); + +int cher_(char *uplo, int *n, float *alpha, fcomplex *x, int *incx, + fcomplex *a, int *lda); + +int cher2_(char *uplo, int *n, fcomplex *alpha, fcomplex *x, int *incx, + fcomplex *y, int *incy, fcomplex *a, int *lda); + +int cher2k_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, + fcomplex *a, int *lda, fcomplex *b, int *ldb, float *beta, + fcomplex *c, int *ldc); + +int cherk_(char *uplo, char *trans, int *n, int *k, float *alpha, + fcomplex *a, int *lda, float *beta, fcomplex *c, int *ldc); + +int chpmv_(char *uplo, int *n, fcomplex *alpha, fcomplex *ap, fcomplex *x, + int *incx, fcomplex *beta, fcomplex *y, int *incy); + +int chpr_(char *uplo, int *n, float *alpha, fcomplex *x, int *incx, + fcomplex *ap); + +int chpr2_(char *uplo, int *n, fcomplex *alpha, fcomplex *x, int *incx, + fcomplex *y, int *incy, fcomplex *ap); + +int crotg_(fcomplex *ca, fcomplex *cb, float *c, fcomplex *s); + +int cscal_(int *n, fcomplex *ca, fcomplex *cx, int *incx); + +int csscal_(int *n, float *sa, fcomplex *cx, int *incx); + +int cswap_(int *n, fcomplex *cx, int *incx, fcomplex *cy, int *incy); + +int csymm_(char *side, char *uplo, int *m, int *n, fcomplex *alpha, + fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, + fcomplex *c, int *ldc); + +int csyr2k_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, + fcomplex *a, int *lda, fcomplex *b, int *ldb, fcomplex *beta, + fcomplex *c, int *ldc); + +int csyrk_(char *uplo, char *trans, int *n, int *k, fcomplex *alpha, + fcomplex *a, int *lda, fcomplex *beta, fcomplex *c, int *ldc); + +int ctbmv_(char *uplo, char *trans, char *diag, int *n, int *k, + fcomplex *a, int *lda, fcomplex *x, int *incx); + +int ctbsv_(char *uplo, char *trans, char *diag, int *n, int *k, + fcomplex *a, int *lda, fcomplex *x, int *incx); + +int ctpmv_(char *uplo, char *trans, char *diag, int *n, fcomplex *ap, + fcomplex *x, int *incx); + +int ctpsv_(char *uplo, char *trans, char *diag, int *n, fcomplex *ap, + fcomplex *x, int *incx); + +int ctrmm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, + int *ldb); + +int ctrmv_(char *uplo, char *trans, char *diag, int *n, fcomplex *a, + int *lda, fcomplex *x, int *incx); + +int ctrsm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, fcomplex *alpha, fcomplex *a, int *lda, fcomplex *b, + int *ldb); + +int ctrsv_(char *uplo, char *trans, char *diag, int *n, fcomplex *a, + int *lda, fcomplex *x, int *incx); + +int daxpy_(int *n, double *sa, double *sx, int *incx, double *sy, + int *incy); + +int dcopy_(int *n, double *sx, int *incx, double *sy, int *incy); + +int dgbmv_(char *trans, int *m, int *n, int *kl, int *ku, + double *alpha, double *a, int *lda, double *x, int *incx, + double *beta, double *y, int *incy); + +int dgemm_(char *transa, char *transb, int *m, int *n, int *k, + double *alpha, double *a, int *lda, double *b, int *ldb, + double *beta, double *c, int *ldc); + +int dgemv_(char *trans, int *m, int *n, double *alpha, double *a, + int *lda, double *x, int *incx, double *beta, double *y, + int *incy); + +int dger_(int *m, int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *a, int *lda); + +int drot_(int *n, double *sx, int *incx, double *sy, int *incy, + double *c, double *s); + +int drotg_(double *sa, double *sb, double *c, double *s); + +int dsbmv_(char *uplo, int *n, int *k, double *alpha, double *a, + int *lda, double *x, int *incx, double *beta, double *y, + int *incy); + +int dscal_(int *n, double *sa, double *sx, int *incx); + +int dspmv_(char *uplo, int *n, double *alpha, double *ap, double *x, + int *incx, double *beta, double *y, int *incy); + +int dspr_(char *uplo, int *n, double *alpha, double *x, int *incx, + double *ap); + +int dspr2_(char *uplo, int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *ap); + +int dswap_(int *n, double *sx, int *incx, double *sy, int *incy); + +int dsymm_(char *side, char *uplo, int *m, int *n, double *alpha, + double *a, int *lda, double *b, int *ldb, double *beta, + double *c, int *ldc); + +int dsymv_(char *uplo, int *n, double *alpha, double *a, int *lda, + double *x, int *incx, double *beta, double *y, int *incy); + +int dsyr_(char *uplo, int *n, double *alpha, double *x, int *incx, + double *a, int *lda); + +int dsyr2_(char *uplo, int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *a, int *lda); + +int dsyr2k_(char *uplo, char *trans, int *n, int *k, double *alpha, + double *a, int *lda, double *b, int *ldb, double *beta, + double *c, int *ldc); + +int dsyrk_(char *uplo, char *trans, int *n, int *k, double *alpha, + double *a, int *lda, double *beta, double *c, int *ldc); + +int dtbmv_(char *uplo, char *trans, char *diag, int *n, int *k, + double *a, int *lda, double *x, int *incx); + +int dtbsv_(char *uplo, char *trans, char *diag, int *n, int *k, + double *a, int *lda, double *x, int *incx); + +int dtpmv_(char *uplo, char *trans, char *diag, int *n, double *ap, + double *x, int *incx); + +int dtpsv_(char *uplo, char *trans, char *diag, int *n, double *ap, + double *x, int *incx); + +int dtrmm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, double *alpha, double *a, int *lda, double *b, + int *ldb); + +int dtrmv_(char *uplo, char *trans, char *diag, int *n, double *a, + int *lda, double *x, int *incx); + +int dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, double *alpha, double *a, int *lda, double *b, + int *ldb); + +int dtrsv_(char *uplo, char *trans, char *diag, int *n, double *a, + int *lda, double *x, int *incx); + + +int saxpy_(int *n, float *sa, float *sx, int *incx, float *sy, int *incy); + +int scopy_(int *n, float *sx, int *incx, float *sy, int *incy); + +int sgbmv_(char *trans, int *m, int *n, int *kl, int *ku, + float *alpha, float *a, int *lda, float *x, int *incx, + float *beta, float *y, int *incy); + +int sgemm_(char *transa, char *transb, int *m, int *n, int *k, + float *alpha, float *a, int *lda, float *b, int *ldb, + float *beta, float *c, int *ldc); + +int sgemv_(char *trans, int *m, int *n, float *alpha, float *a, + int *lda, float *x, int *incx, float *beta, float *y, + int *incy); + +int sger_(int *m, int *n, float *alpha, float *x, int *incx, + float *y, int *incy, float *a, int *lda); + +int srot_(int *n, float *sx, int *incx, float *sy, int *incy, + float *c, float *s); + +int srotg_(float *sa, float *sb, float *c, float *s); + +int ssbmv_(char *uplo, int *n, int *k, float *alpha, float *a, + int *lda, float *x, int *incx, float *beta, float *y, + int *incy); + +int sscal_(int *n, float *sa, float *sx, int *incx); + +int sspmv_(char *uplo, int *n, float *alpha, float *ap, float *x, + int *incx, float *beta, float *y, int *incy); + +int sspr_(char *uplo, int *n, float *alpha, float *x, int *incx, + float *ap); + +int sspr2_(char *uplo, int *n, float *alpha, float *x, int *incx, + float *y, int *incy, float *ap); + +int sswap_(int *n, float *sx, int *incx, float *sy, int *incy); + +int ssymm_(char *side, char *uplo, int *m, int *n, float *alpha, + float *a, int *lda, float *b, int *ldb, float *beta, + float *c, int *ldc); + +int ssymv_(char *uplo, int *n, float *alpha, float *a, int *lda, + float *x, int *incx, float *beta, float *y, int *incy); + +int ssyr_(char *uplo, int *n, float *alpha, float *x, int *incx, + float *a, int *lda); + +int ssyr2_(char *uplo, int *n, float *alpha, float *x, int *incx, + float *y, int *incy, float *a, int *lda); + +int ssyr2k_(char *uplo, char *trans, int *n, int *k, float *alpha, + float *a, int *lda, float *b, int *ldb, float *beta, + float *c, int *ldc); + +int ssyrk_(char *uplo, char *trans, int *n, int *k, float *alpha, + float *a, int *lda, float *beta, float *c, int *ldc); + +int stbmv_(char *uplo, char *trans, char *diag, int *n, int *k, + float *a, int *lda, float *x, int *incx); + +int stbsv_(char *uplo, char *trans, char *diag, int *n, int *k, + float *a, int *lda, float *x, int *incx); + +int stpmv_(char *uplo, char *trans, char *diag, int *n, float *ap, + float *x, int *incx); + +int stpsv_(char *uplo, char *trans, char *diag, int *n, float *ap, + float *x, int *incx); + +int strmm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, float *alpha, float *a, int *lda, float *b, + int *ldb); + +int strmv_(char *uplo, char *trans, char *diag, int *n, float *a, + int *lda, float *x, int *incx); + +int strsm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, float *alpha, float *a, int *lda, float *b, + int *ldb); + +int strsv_(char *uplo, char *trans, char *diag, int *n, float *a, + int *lda, float *x, int *incx); + +int zaxpy_(int *n, dcomplex *ca, dcomplex *cx, int *incx, dcomplex *cy, + int *incy); + +int zcopy_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); + +int zdscal_(int *n, double *sa, dcomplex *cx, int *incx); + +int zgbmv_(char *trans, int *m, int *n, int *kl, int *ku, + dcomplex *alpha, dcomplex *a, int *lda, dcomplex *x, int *incx, + dcomplex *beta, dcomplex *y, int *incy); + +int zgemm_(char *transa, char *transb, int *m, int *n, int *k, + dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, int *ldb, + dcomplex *beta, dcomplex *c, int *ldc); + +int zgemv_(char *trans, int *m, int *n, dcomplex *alpha, dcomplex *a, + int *lda, dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, + int *incy); + +int zgerc_(int *m, int *n, dcomplex *alpha, dcomplex *x, int *incx, + dcomplex *y, int *incy, dcomplex *a, int *lda); + +int zgeru_(int *m, int *n, dcomplex *alpha, dcomplex *x, int *incx, + dcomplex *y, int *incy, dcomplex *a, int *lda); + +int zhbmv_(char *uplo, int *n, int *k, dcomplex *alpha, dcomplex *a, + int *lda, dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, + int *incy); + +int zhemm_(char *side, char *uplo, int *m, int *n, dcomplex *alpha, + dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, + dcomplex *c, int *ldc); + +int zhemv_(char *uplo, int *n, dcomplex *alpha, dcomplex *a, int *lda, + dcomplex *x, int *incx, dcomplex *beta, dcomplex *y, int *incy); + +int zher_(char *uplo, int *n, double *alpha, dcomplex *x, int *incx, + dcomplex *a, int *lda); + +int zher2_(char *uplo, int *n, dcomplex *alpha, dcomplex *x, int *incx, + dcomplex *y, int *incy, dcomplex *a, int *lda); + +int zher2k_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, + dcomplex *a, int *lda, dcomplex *b, int *ldb, double *beta, + dcomplex *c, int *ldc); + +int zherk_(char *uplo, char *trans, int *n, int *k, double *alpha, + dcomplex *a, int *lda, double *beta, dcomplex *c, int *ldc); + +int zhpmv_(char *uplo, int *n, dcomplex *alpha, dcomplex *ap, dcomplex *x, + int *incx, dcomplex *beta, dcomplex *y, int *incy); + +int zhpr_(char *uplo, int *n, double *alpha, dcomplex *x, int *incx, + dcomplex *ap); + +int zhpr2_(char *uplo, int *n, dcomplex *alpha, dcomplex *x, int *incx, + dcomplex *y, int *incy, dcomplex *ap); + +int zrotg_(dcomplex *ca, dcomplex *cb, double *c, dcomplex *s); + +int zscal_(int *n, dcomplex *ca, dcomplex *cx, int *incx); + +int zswap_(int *n, dcomplex *cx, int *incx, dcomplex *cy, int *incy); + +int zsymm_(char *side, char *uplo, int *m, int *n, dcomplex *alpha, + dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, + dcomplex *c, int *ldc); + +int zsyr2k_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, + dcomplex *a, int *lda, dcomplex *b, int *ldb, dcomplex *beta, + dcomplex *c, int *ldc); + +int zsyrk_(char *uplo, char *trans, int *n, int *k, dcomplex *alpha, + dcomplex *a, int *lda, dcomplex *beta, dcomplex *c, int *ldc); + +int ztbmv_(char *uplo, char *trans, char *diag, int *n, int *k, + dcomplex *a, int *lda, dcomplex *x, int *incx); + +int ztbsv_(char *uplo, char *trans, char *diag, int *n, int *k, + dcomplex *a, int *lda, dcomplex *x, int *incx); + +int ztpmv_(char *uplo, char *trans, char *diag, int *n, dcomplex *ap, + dcomplex *x, int *incx); + +int ztpsv_(char *uplo, char *trans, char *diag, int *n, dcomplex *ap, + dcomplex *x, int *incx); + +int ztrmm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, + int *ldb); + +int ztrmv_(char *uplo, char *trans, char *diag, int *n, dcomplex *a, + int *lda, dcomplex *x, int *incx); + +int ztrsm_(char *side, char *uplo, char *transa, char *diag, int *m, + int *n, dcomplex *alpha, dcomplex *a, int *lda, dcomplex *b, + int *ldb); + +int ztrsv_(char *uplo, char *trans, char *diag, int *n, dcomplex *a, + int *lda, dcomplex *x, int *incx); diff --git a/random/include/LibLinear/blas/daxpy.c b/random/include/LibLinear/blas/daxpy.c new file mode 100755 index 0000000..62c50bd --- /dev/null +++ b/random/include/LibLinear/blas/daxpy.c @@ -0,0 +1,49 @@ +#include "blas.h" + +int daxpy_(int *n, double *sa, double *sx, int *incx, double *sy, + int *incy) +{ + long int i, m, ix, iy, nn, iincx, iincy; + register double ssa; + + /* constant times a vector plus a vector. + uses unrolled loop for increments equal to one. + jack dongarra, linpack, 3/11/78. + modified 12/3/93, array(1) declarations changed to array(*) */ + + /* Dereference inputs */ + nn = *n; + ssa = *sa; + iincx = *incx; + iincy = *incy; + + if( nn > 0 && ssa != 0.0 ) + { + if (iincx == 1 && iincy == 1) /* code for both increments equal to 1 */ + { + m = nn-3; + for (i = 0; i < m; i += 4) + { + sy[i] += ssa * sx[i]; + sy[i+1] += ssa * sx[i+1]; + sy[i+2] += ssa * sx[i+2]; + sy[i+3] += ssa * sx[i+3]; + } + for ( ; i < nn; ++i) /* clean-up loop */ + sy[i] += ssa * sx[i]; + } + else /* code for unequal increments or equal increments not equal to 1 */ + { + ix = iincx >= 0 ? 0 : (1 - nn) * iincx; + iy = iincy >= 0 ? 0 : (1 - nn) * iincy; + for (i = 0; i < nn; i++) + { + sy[iy] += ssa * sx[ix]; + ix += iincx; + iy += iincy; + } + } + } + + return 0; +} /* daxpy_ */ diff --git a/random/include/LibLinear/blas/ddot.c b/random/include/LibLinear/blas/ddot.c new file mode 100755 index 0000000..e77c7c1 --- /dev/null +++ b/random/include/LibLinear/blas/ddot.c @@ -0,0 +1,50 @@ +#include "blas.h" + +double ddot_(int *n, double *sx, int *incx, double *sy, int *incy) +{ + long int i, m, nn, iincx, iincy; + double stemp; + long int ix, iy; + + /* forms the dot product of two vectors. + uses unrolled loops for increments equal to one. + jack dongarra, linpack, 3/11/78. + modified 12/3/93, array(1) declarations changed to array(*) */ + + /* Dereference inputs */ + nn = *n; + iincx = *incx; + iincy = *incy; + + stemp = 0.0; + if (nn > 0) + { + if (iincx == 1 && iincy == 1) /* code for both increments equal to 1 */ + { + m = nn-4; + for (i = 0; i < m; i += 5) + stemp += sx[i] * sy[i] + sx[i+1] * sy[i+1] + sx[i+2] * sy[i+2] + + sx[i+3] * sy[i+3] + sx[i+4] * sy[i+4]; + + for ( ; i < nn; i++) /* clean-up loop */ + stemp += sx[i] * sy[i]; + } + else /* code for unequal increments or equal increments not equal to 1 */ + { + ix = 0; + iy = 0; + if (iincx < 0) + ix = (1 - nn) * iincx; + if (iincy < 0) + iy = (1 - nn) * iincy; + for (i = 0; i < nn; i++) + { + stemp += sx[ix] * sy[iy]; + ix += iincx; + iy += iincy; + } + } + } + + return stemp; +} /* ddot_ */ diff --git a/random/include/LibLinear/blas/dnrm2.c b/random/include/LibLinear/blas/dnrm2.c new file mode 100755 index 0000000..2facb8e --- /dev/null +++ b/random/include/LibLinear/blas/dnrm2.c @@ -0,0 +1,62 @@ +#include /* Needed for fabs() and sqrt() */ +#include "blas.h" + +double dnrm2_(int *n, double *x, int *incx) +{ + long int ix, nn, iincx; + double norm, scale, absxi, ssq, temp; + +/* DNRM2 returns the euclidean norm of a vector via the function + name, so that + + DNRM2 := sqrt( x'*x ) + + -- This version written on 25-October-1982. + Modified on 14-October-1993 to inline the call to SLASSQ. + Sven Hammarling, Nag Ltd. */ + + /* Dereference inputs */ + nn = *n; + iincx = *incx; + + if( nn > 0 && iincx > 0 ) + { + if (nn == 1) + { + norm = fabs(x[0]); + } + else + { + scale = 0.0; + ssq = 1.0; + + /* The following loop is equivalent to this call to the LAPACK + auxiliary routine: CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */ + + for (ix=(nn-1)*iincx; ix>=0; ix-=iincx) + { + if (x[ix] != 0.0) + { + absxi = fabs(x[ix]); + if (scale < absxi) + { + temp = scale / absxi; + ssq = ssq * (temp * temp) + 1.0; + scale = absxi; + } + else + { + temp = absxi / scale; + ssq += temp * temp; + } + } + } + norm = scale * sqrt(ssq); + } + } + else + norm = 0.0; + + return norm; + +} /* dnrm2_ */ diff --git a/random/include/LibLinear/blas/dscal.c b/random/include/LibLinear/blas/dscal.c new file mode 100755 index 0000000..ee73921 --- /dev/null +++ b/random/include/LibLinear/blas/dscal.c @@ -0,0 +1,44 @@ +#include "blas.h" + +int dscal_(int *n, double *sa, double *sx, int *incx) +{ + long int i, m, nincx, nn, iincx; + double ssa; + + /* scales a vector by a constant. + uses unrolled loops for increment equal to 1. + jack dongarra, linpack, 3/11/78. + modified 3/93 to return if incx .le. 0. + modified 12/3/93, array(1) declarations changed to array(*) */ + + /* Dereference inputs */ + nn = *n; + iincx = *incx; + ssa = *sa; + + if (nn > 0 && iincx > 0) + { + if (iincx == 1) /* code for increment equal to 1 */ + { + m = nn-4; + for (i = 0; i < m; i += 5) + { + sx[i] = ssa * sx[i]; + sx[i+1] = ssa * sx[i+1]; + sx[i+2] = ssa * sx[i+2]; + sx[i+3] = ssa * sx[i+3]; + sx[i+4] = ssa * sx[i+4]; + } + for ( ; i < nn; ++i) /* clean-up loop */ + sx[i] = ssa * sx[i]; + } + else /* code for increment not equal to 1 */ + { + nincx = nn * iincx; + for (i = 0; i < nincx; i += iincx) + sx[i] = ssa * sx[i]; + } + } + + return 0; +} /* dscal_ */ diff --git a/random/include/LibLinear/linear.cpp b/random/include/LibLinear/linear.cpp new file mode 100755 index 0000000..837a565 --- /dev/null +++ b/random/include/LibLinear/linear.cpp @@ -0,0 +1,2811 @@ +#include +#include +#include +#include +#include +#include +#include "linear.h" +#include "tron.h" +typedef signed char schar; +template static inline void swap(T& x, T& y) { T t=x; x=y; y=t; } +#ifndef min +template static inline T min(T x,T y) { return (x static inline T max(T x,T y) { return (x>y)?x:y; } +#endif +template static inline void clone(T*& dst, S* src, int n) +{ + dst = new T[n]; + memcpy((void *)dst,(void *)src,sizeof(T)*n); +} +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) +#define INF HUGE_VAL + +static void print_string_stdout(const char *s) +{ + fputs(s,stdout); + fflush(stdout); +} + +static void (*liblinear_print_string) (const char *) = &print_string_stdout; + +#if 1 +static void info(const char *fmt,...) +{ + char buf[BUFSIZ]; + va_list ap; + va_start(ap,fmt); + vsprintf(buf,fmt,ap); + va_end(ap); + (*liblinear_print_string)(buf); +} +#else +static void info(const char *fmt,...) {} +#endif + +class l2r_lr_fun: public function +{ +public: + l2r_lr_fun(const problem *prob, double *C); + ~l2r_lr_fun(); + + double fun(double *w); + void grad(double *w, double *g); + void Hv(double *s, double *Hs); + + int get_nr_variable(void); + +private: + void Xv(double *v, double *Xv); + void XTv(double *v, double *XTv); + + double *C; + double *z; + double *D; + const problem *prob; +}; + +l2r_lr_fun::l2r_lr_fun(const problem *prob, double *C) +{ + int l=prob->l; + + this->prob = prob; + + z = new double[l]; + D = new double[l]; + this->C = C; +} + +l2r_lr_fun::~l2r_lr_fun() +{ + delete[] z; + delete[] D; +} + + +double l2r_lr_fun::fun(double *w) +{ + int i; + double f=0; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + + Xv(w, z); + + for(i=0;i= 0) + f += C[i]*log(1 + exp(-yz)); + else + f += C[i]*(-yz+log(1 + exp(yz))); + } + + return(f); +} + +void l2r_lr_fun::grad(double *w, double *g) +{ + int i; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + + for(i=0;in; +} + +void l2r_lr_fun::Hv(double *s, double *Hs) +{ + int i; + int l=prob->l; + int w_size=get_nr_variable(); + double *wa = new double[l]; + + Xv(s, wa); + for(i=0;il; + feature_node **x=prob->x; + + for(i=0;iindex!=-1) + { + Xv[i]+=v[s->index-1]*s->value; + s++; + } + } +} + +void l2r_lr_fun::XTv(double *v, double *XTv) +{ + int i; + int l=prob->l; + int w_size=get_nr_variable(); + feature_node **x=prob->x; + + for(i=0;iindex!=-1) + { + XTv[s->index-1]+=v[i]*s->value; + s++; + } + } +} + +class l2r_l2_svc_fun: public function +{ +public: + l2r_l2_svc_fun(const problem *prob, double *C); + ~l2r_l2_svc_fun(); + + double fun(double *w); + void grad(double *w, double *g); + void Hv(double *s, double *Hs); + + int get_nr_variable(void); + +protected: + void Xv(double *v, double *Xv); + void subXv(double *v, double *Xv); + void subXTv(double *v, double *XTv); + + double *C; + double *z; + double *D; + int *I; + int sizeI; + const problem *prob; +}; + +l2r_l2_svc_fun::l2r_l2_svc_fun(const problem *prob, double *C) +{ + int l=prob->l; + + this->prob = prob; + + z = new double[l]; + D = new double[l]; + I = new int[l]; + this->C = C; +} + +l2r_l2_svc_fun::~l2r_l2_svc_fun() +{ + delete[] z; + delete[] D; + delete[] I; +} + +double l2r_l2_svc_fun::fun(double *w) +{ + int i; + double f=0; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + + Xv(w, z); + + for(i=0;i 0) + f += C[i]*d*d; + } + + return(f); +} + +void l2r_l2_svc_fun::grad(double *w, double *g) +{ + int i; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + + sizeI = 0; + for (i=0;in; +} + +void l2r_l2_svc_fun::Hv(double *s, double *Hs) +{ + int i; + int w_size=get_nr_variable(); + double *wa = new double[sizeI]; + + subXv(s, wa); + for(i=0;il; + feature_node **x=prob->x; + + for(i=0;iindex!=-1) + { + Xv[i]+=v[s->index-1]*s->value; + s++; + } + } +} + +void l2r_l2_svc_fun::subXv(double *v, double *Xv) +{ + int i; + feature_node **x=prob->x; + + for(i=0;iindex!=-1) + { + Xv[i]+=v[s->index-1]*s->value; + s++; + } + } +} + +void l2r_l2_svc_fun::subXTv(double *v, double *XTv) +{ + int i; + int w_size=get_nr_variable(); + feature_node **x=prob->x; + + for(i=0;iindex!=-1) + { + XTv[s->index-1]+=v[i]*s->value; + s++; + } + } +} + +class l2r_l2_svr_fun: public l2r_l2_svc_fun +{ +public: + l2r_l2_svr_fun(const problem *prob, double *C, double p); + + double fun(double *w); + void grad(double *w, double *g); + +private: + double p; +}; + +l2r_l2_svr_fun::l2r_l2_svr_fun(const problem *prob, double *C, double p): + l2r_l2_svc_fun(prob, C) +{ + this->p = p; +} + +double l2r_l2_svr_fun::fun(double *w) +{ + int i; + double f=0; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + double d; + + Xv(w, z); + + for(i=0;i p) + f += C[i]*(d-p)*(d-p); + } + + return(f); +} + +void l2r_l2_svr_fun::grad(double *w, double *g) +{ + int i; + double *y=prob->y; + int l=prob->l; + int w_size=get_nr_variable(); + double d; + + sizeI = 0; + for(i=0;i p) + { + z[sizeI] = C[i]*(d-p); + I[sizeI] = i; + sizeI++; + } + + } + subXTv(z, g); + + for(i=0;iy[i]) +// To support weights for instances, use GETI(i) (i) + +class Solver_MCSVM_CS +{ + public: + Solver_MCSVM_CS(const problem *prob, int nr_class, double *C, double eps=0.1, int max_iter=100000); + ~Solver_MCSVM_CS(); + void Solve(double *w); + private: + void solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new); + bool be_shrunk(int i, int m, int yi, double alpha_i, double minG); + double *B, *C, *G; + int w_size, l; + int nr_class; + int max_iter; + double eps; + const problem *prob; +}; + +Solver_MCSVM_CS::Solver_MCSVM_CS(const problem *prob, int nr_class, double *weighted_C, double eps, int max_iter) +{ + this->w_size = prob->n; + this->l = prob->l; + this->nr_class = nr_class; + this->eps = eps; + this->max_iter = max_iter; + this->prob = prob; + this->B = new double[nr_class]; + this->G = new double[nr_class]; + this->C = weighted_C; +} + +Solver_MCSVM_CS::~Solver_MCSVM_CS() +{ + delete[] B; + delete[] G; +} + +int compare_double(const void *a, const void *b) +{ + if(*(double *)a > *(double *)b) + return -1; + if(*(double *)a < *(double *)b) + return 1; + return 0; +} + +void Solver_MCSVM_CS::solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new) +{ + int r; + double *D; + + clone(D, B, active_i); + if(yi < active_i) + D[yi] += A_i*C_yi; + qsort(D, active_i, sizeof(double), compare_double); + + double beta = D[0] - A_i*C_yi; + for(r=1;ry[i] == m + // alpha[i*nr_class+m] <= 0 if prob->y[i] != m + // If initial alpha isn't zero, uncomment the for loop below to initialize w + for(i=0;ix[i]; + QD[i] = 0; + while(xi->index != -1) + { + double val = xi->value; + QD[i] += val*val; + + // Uncomment the for loop if initial alpha isn't zero + // for(m=0; mindex-1)*nr_class+m] += alpha[i*nr_class+m]*val; + xi++; + } + active_size_i[i] = nr_class; + y_index[i] = (int)prob->y[i]; + index[i] = i; + } + + while(iter < max_iter) + { + double stopping = -INF; + for(i=0;i 0) + { + for(m=0;mx[i]; + while(xi->index!= -1) + { + double *w_i = &w[(xi->index-1)*nr_class]; + for(m=0;mvalue); + xi++; + } + + double minG = INF; + double maxG = -INF; + for(m=0;m maxG) + maxG = G[m]; + } + if(y_index[i] < active_size_i[i]) + if(alpha_i[(int) prob->y[i]] < C[GETI(i)] && G[y_index[i]] < minG) + minG = G[y_index[i]]; + + for(m=0;mm) + { + if(!be_shrunk(i, active_size_i[i], y_index[i], + alpha_i[alpha_index_i[active_size_i[i]]], minG)) + { + swap(alpha_index_i[m], alpha_index_i[active_size_i[i]]); + swap(G[m], G[active_size_i[i]]); + if(y_index[i] == active_size_i[i]) + y_index[i] = m; + else if(y_index[i] == m) + y_index[i] = active_size_i[i]; + break; + } + active_size_i[i]--; + } + } + } + + if(active_size_i[i] <= 1) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + + if(maxG-minG <= 1e-12) + continue; + else + stopping = max(maxG - minG, stopping); + + for(m=0;m= 1e-12) + { + d_ind[nz_d] = alpha_index_i[m]; + d_val[nz_d] = d; + nz_d++; + } + } + + xi = prob->x[i]; + while(xi->index != -1) + { + double *w_i = &w[(xi->index-1)*nr_class]; + for(m=0;mvalue; + xi++; + } + } + } + + iter++; + if(iter % 10 == 0) + { + info("."); + } + + if(stopping < eps_shrink) + { + if(stopping < eps && start_from_all == true) + break; + else + { + active_size = l; + for(i=0;i= max_iter) + info("\nWARNING: reaching max number of iterations\n"); + + // calculate objective value + double v = 0; + int nSV = 0; + for(i=0;i 0) + nSV++; + } + for(i=0;iy[i]]; + info("Objective value = %lf\n",v); + info("nSV = %d\n",nSV); + + delete [] alpha; + delete [] alpha_new; + delete [] index; + delete [] QD; + delete [] d_ind; + delete [] d_val; + delete [] alpha_index; + delete [] y_index; + delete [] active_size_i; +} + +// A coordinate descent algorithm for +// L1-loss and L2-loss SVM dual problems +// +// min_\alpha 0.5(\alpha^T (Q + D)\alpha) - e^T \alpha, +// s.t. 0 <= \alpha_i <= upper_bound_i, +// +// where Qij = yi yj xi^T xj and +// D is a diagonal matrix +// +// In L1-SVM case: +// upper_bound_i = Cp if y_i = 1 +// upper_bound_i = Cn if y_i = -1 +// D_ii = 0 +// In L2-SVM case: +// upper_bound_i = INF +// D_ii = 1/(2*Cp) if y_i = 1 +// D_ii = 1/(2*Cn) if y_i = -1 +// +// Given: +// x, y, Cp, Cn +// eps is the stopping tolerance +// +// solution will be put in w +// +// See Algorithm 3 of Hsieh et al., ICML 2008 + +#undef GETI +#define GETI(i) (y[i]+1) +// To support weights for instances, use GETI(i) (i) + +static void solve_l2r_l1l2_svc( + const problem *prob, double *w, double eps, + double Cp, double Cn, int solver_type) +{ + int l = prob->l; + int w_size = prob->n; + int i, s, iter = 0; + double C, d, G; + double *QD = new double[l]; + int max_iter = 1000; + int *index = new int[l]; + double *alpha = new double[l]; + schar *y = new schar[l]; + int active_size = l; + + // PG: projected gradient, for shrinking and stopping + double PG; + double PGmax_old = INF; + double PGmin_old = -INF; + double PGmax_new, PGmin_new; + + // default solver_type: L2R_L2LOSS_SVC_DUAL + double diag[3] = {0.5/Cn, 0, 0.5/Cp}; + double upper_bound[3] = {INF, 0, INF}; + if(solver_type == L2R_L1LOSS_SVC_DUAL) + { + diag[0] = 0; + diag[2] = 0; + upper_bound[0] = Cn; + upper_bound[2] = Cp; + } + + for(i=0; iy[i] > 0) + { + y[i] = +1; + } + else + { + y[i] = -1; + } + } + + // Initial alpha can be set here. Note that + // 0 <= alpha[i] <= upper_bound[GETI(i)] + for(i=0; ix[i]; + while (xi->index != -1) + { + double val = xi->value; + QD[i] += val*val; + w[xi->index-1] += y[i]*alpha[i]*val; + xi++; + } + index[i] = i; + } + + while (iter < max_iter) + { + PGmax_new = -INF; + PGmin_new = INF; + + for (i=0; ix[i]; + while(xi->index!= -1) + { + G += w[xi->index-1]*(xi->value); + xi++; + } + G = G*yi-1; + + C = upper_bound[GETI(i)]; + G += alpha[i]*diag[GETI(i)]; + + PG = 0; + if (alpha[i] == 0) + { + if (G > PGmax_old) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + else if (G < 0) + PG = G; + } + else if (alpha[i] == C) + { + if (G < PGmin_old) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + else if (G > 0) + PG = G; + } + else + PG = G; + + PGmax_new = max(PGmax_new, PG); + PGmin_new = min(PGmin_new, PG); + + if(fabs(PG) > 1.0e-12) + { + double alpha_old = alpha[i]; + alpha[i] = min(max(alpha[i] - G/QD[i], 0.0), C); + d = (alpha[i] - alpha_old)*yi; + xi = prob->x[i]; + while (xi->index != -1) + { + w[xi->index-1] += d*xi->value; + xi++; + } + } + } + + iter++; + if(iter % 10 == 0) + info("."); + + if(PGmax_new - PGmin_new <= eps) + { + if(active_size == l) + break; + else + { + active_size = l; + info("*"); + PGmax_old = INF; + PGmin_old = -INF; + continue; + } + } + PGmax_old = PGmax_new; + PGmin_old = PGmin_new; + if (PGmax_old <= 0) + PGmax_old = INF; + if (PGmin_old >= 0) + PGmin_old = -INF; + } + + info("\noptimization finished, #iter = %d\n",iter); + if (iter >= max_iter) + info("\nWARNING: reaching max number of iterations\nUsing -s 2 may be faster (also see FAQ)\n\n"); + + // calculate objective value + + double v = 0; + int nSV = 0; + for(i=0; i 0) + ++nSV; + } + info("Objective value = %lf\n",v/2); + info("nSV = %d\n",nSV); + + delete [] QD; + delete [] alpha; + delete [] y; + delete [] index; +} + + +// A coordinate descent algorithm for +// L1-loss and L2-loss epsilon-SVR dual problem +// +// min_\beta 0.5\beta^T (Q + diag(lambda)) \beta - p \sum_{i=1}^l|\beta_i| + \sum_{i=1}^l yi\beta_i, +// s.t. -upper_bound_i <= \beta_i <= upper_bound_i, +// +// where Qij = xi^T xj and +// D is a diagonal matrix +// +// In L1-SVM case: +// upper_bound_i = C +// lambda_i = 0 +// In L2-SVM case: +// upper_bound_i = INF +// lambda_i = 1/(2*C) +// +// Given: +// x, y, p, C +// eps is the stopping tolerance +// +// solution will be put in w +// +// See Algorithm 4 of Ho and Lin, 2012 + +#undef GETI +#define GETI(i) (0) +// To support weights for instances, use GETI(i) (i) + +static void solve_l2r_l1l2_svr( + const problem *prob, double *w, const parameter *param, + int solver_type) +{ + int l = prob->l; + double C = param->C; + double p = param->p; + int w_size = prob->n; + double eps = param->eps; + int i, s, iter = 0; + int max_iter = 1000; + int active_size = l; + int *index = new int[l]; + + double d, G, H; + double Gmax_old = INF; + double Gmax_new, Gnorm1_new; + double Gnorm1_init; + double *beta = new double[l]; + double *QD = new double[l]; + double *y = prob->y; + + // L2R_L2LOSS_SVR_DUAL + double lambda[1], upper_bound[1]; + lambda[0] = 0.5/C; + upper_bound[0] = INF; + + if(solver_type == L2R_L1LOSS_SVR_DUAL) + { + lambda[0] = 0; + upper_bound[0] = C; + } + + // Initial beta can be set here. Note that + // -upper_bound <= beta[i] <= upper_bound + for(i=0; ix[i]; + while(xi->index != -1) + { + double val = xi->value; + QD[i] += val*val; + w[xi->index-1] += beta[i]*val; + xi++; + } + + index[i] = i; + } + + + while(iter < max_iter) + { + Gmax_new = 0; + Gnorm1_new = 0; + + for(i=0; ix[i]; + while(xi->index != -1) + { + int ind = xi->index-1; + double val = xi->value; + G += val*w[ind]; + xi++; + } + + double Gp = G+p; + double Gn = G-p; + double violation = 0; + if(beta[i] == 0) + { + if(Gp < 0) + violation = -Gp; + else if(Gn > 0) + violation = Gn; + else if(Gp>Gmax_old && Gn<-Gmax_old) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + } + else if(beta[i] >= upper_bound[GETI(i)]) + { + if(Gp > 0) + violation = Gp; + else if(Gp < -Gmax_old) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + } + else if(beta[i] <= -upper_bound[GETI(i)]) + { + if(Gn < 0) + violation = -Gn; + else if(Gn > Gmax_old) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + } + else if(beta[i] > 0) + violation = fabs(Gp); + else + violation = fabs(Gn); + + Gmax_new = max(Gmax_new, violation); + Gnorm1_new += violation; + + // obtain Newton direction d + if(Gp < H*beta[i]) + d = -Gp/H; + else if(Gn > H*beta[i]) + d = -Gn/H; + else + d = -beta[i]; + + if(fabs(d) < 1.0e-12) + continue; + + double beta_old = beta[i]; + beta[i] = min(max(beta[i]+d, -upper_bound[GETI(i)]), upper_bound[GETI(i)]); + d = beta[i]-beta_old; + + if(d != 0) + { + xi = prob->x[i]; + while(xi->index != -1) + { + w[xi->index-1] += d*xi->value; + xi++; + } + } + } + + if(iter == 0) + Gnorm1_init = Gnorm1_new; + iter++; + if(iter % 10 == 0) + info("."); + + if(Gnorm1_new <= eps*Gnorm1_init) + { + if(active_size == l) + break; + else + { + active_size = l; + info("*"); + Gmax_old = INF; + continue; + } + } + + Gmax_old = Gmax_new; + } + + info("\noptimization finished, #iter = %d\n", iter); + if(iter >= max_iter) + info("\nWARNING: reaching max number of iterations\nUsing -s 11 may be faster\n\n"); + + // calculate objective value + double v = 0; + int nSV = 0; + for(i=0; il; + int w_size = prob->n; + int i, s, iter = 0; + double *xTx = new double[l]; + int max_iter = 1000; + int *index = new int[l]; + double *alpha = new double[2*l]; // store alpha and C - alpha + schar *y = new schar[l]; + int max_inner_iter = 100; // for inner Newton + double innereps = 1e-2; + double innereps_min = min(1e-8, eps); + double upper_bound[3] = {Cn, 0, Cp}; + + for(i=0; iy[i] > 0) + { + y[i] = +1; + } + else + { + y[i] = -1; + } + } + + // Initial alpha can be set here. Note that + // 0 < alpha[i] < upper_bound[GETI(i)] + // alpha[2*i] + alpha[2*i+1] = upper_bound[GETI(i)] + for(i=0; ix[i]; + while (xi->index != -1) + { + double val = xi->value; + xTx[i] += val*val; + w[xi->index-1] += y[i]*alpha[2*i]*val; + xi++; + } + index[i] = i; + } + + while (iter < max_iter) + { + for (i=0; ix[i]; + while (xi->index != -1) + { + ywTx += w[xi->index-1]*xi->value; + xi++; + } + ywTx *= y[i]; + double a = xisq, b = ywTx; + + // Decide to minimize g_1(z) or g_2(z) + int ind1 = 2*i, ind2 = 2*i+1, sign = 1; + if(0.5*a*(alpha[ind2]-alpha[ind1])+b < 0) + { + ind1 = 2*i+1; + ind2 = 2*i; + sign = -1; + } + + // g_t(z) = z*log(z) + (C-z)*log(C-z) + 0.5a(z-alpha_old)^2 + sign*b(z-alpha_old) + double alpha_old = alpha[ind1]; + double z = alpha_old; + if(C - z < 0.5 * C) + z = 0.1*z; + double gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); + Gmax = max(Gmax, fabs(gp)); + + // Newton method on the sub-problem + const double eta = 0.1; // xi in the paper + int inner_iter = 0; + while (inner_iter <= max_inner_iter) + { + if(fabs(gp) < innereps) + break; + double gpp = a + C/(C-z)/z; + double tmpz = z - gp/gpp; + if(tmpz <= 0) + z *= eta; + else // tmpz in (0, C) + z = tmpz; + gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); + newton_iter++; + inner_iter++; + } + + if(inner_iter > 0) // update w + { + alpha[ind1] = z; + alpha[ind2] = C-z; + xi = prob->x[i]; + while (xi->index != -1) + { + w[xi->index-1] += sign*(z-alpha_old)*yi*xi->value; + xi++; + } + } + } + + iter++; + if(iter % 10 == 0) + info("."); + + if(Gmax < eps) + break; + + if(newton_iter <= l/10) + innereps = max(innereps_min, 0.1*innereps); + + } + + info("\noptimization finished, #iter = %d\n",iter); + if (iter >= max_iter) + info("\nWARNING: reaching max number of iterations\nUsing -s 0 may be faster (also see FAQ)\n\n"); + + // calculate objective value + + double v = 0; + for(i=0; il; + int w_size = prob_col->n; + int j, s, iter = 0; + int max_iter = 1000; + int active_size = w_size; + int max_num_linesearch = 20; + + double sigma = 0.01; + double d, G_loss, G, H; + double Gmax_old = INF; + double Gmax_new, Gnorm1_new; + double Gnorm1_init; + double d_old, d_diff; + double loss_old, loss_new; + double appxcond, cond; + + int *index = new int[w_size]; + schar *y = new schar[l]; + double *b = new double[l]; // b = 1-ywTx + double *xj_sq = new double[w_size]; + feature_node *x; + + double C[3] = {Cn,0,Cp}; + + // Initial w can be set here. + for(j=0; jy[j] > 0) + y[j] = 1; + else + y[j] = -1; + } + for(j=0; jx[j]; + while(x->index != -1) + { + int ind = x->index-1; + x->value *= y[ind]; // x->value stores yi*xij + double val = x->value; + b[ind] -= w[j]*val; + xj_sq[j] += C[GETI(ind)]*val*val; + x++; + } + } + + while(iter < max_iter) + { + Gmax_new = 0; + Gnorm1_new = 0; + + for(j=0; jx[j]; + while(x->index != -1) + { + int ind = x->index-1; + if(b[ind] > 0) + { + double val = x->value; + double tmp = C[GETI(ind)]*val; + G_loss -= tmp*b[ind]; + H += tmp*val; + } + x++; + } + G_loss *= 2; + + G = G_loss; + H *= 2; + H = max(H, 1e-12); + + double Gp = G+1; + double Gn = G-1; + double violation = 0; + if(w[j] == 0) + { + if(Gp < 0) + violation = -Gp; + else if(Gn > 0) + violation = Gn; + else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + } + else if(w[j] > 0) + violation = fabs(Gp); + else + violation = fabs(Gn); + + Gmax_new = max(Gmax_new, violation); + Gnorm1_new += violation; + + // obtain Newton direction d + if(Gp < H*w[j]) + d = -Gp/H; + else if(Gn > H*w[j]) + d = -Gn/H; + else + d = -w[j]; + + if(fabs(d) < 1.0e-12) + continue; + + double delta = fabs(w[j]+d)-fabs(w[j]) + G*d; + d_old = 0; + int num_linesearch; + for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++) + { + d_diff = d_old - d; + cond = fabs(w[j]+d)-fabs(w[j]) - sigma*delta; + + appxcond = xj_sq[j]*d*d + G_loss*d + cond; + if(appxcond <= 0) + { + x = prob_col->x[j]; + while(x->index != -1) + { + b[x->index-1] += d_diff*x->value; + x++; + } + break; + } + + if(num_linesearch == 0) + { + loss_old = 0; + loss_new = 0; + x = prob_col->x[j]; + while(x->index != -1) + { + int ind = x->index-1; + if(b[ind] > 0) + loss_old += C[GETI(ind)]*b[ind]*b[ind]; + double b_new = b[ind] + d_diff*x->value; + b[ind] = b_new; + if(b_new > 0) + loss_new += C[GETI(ind)]*b_new*b_new; + x++; + } + } + else + { + loss_new = 0; + x = prob_col->x[j]; + while(x->index != -1) + { + int ind = x->index-1; + double b_new = b[ind] + d_diff*x->value; + b[ind] = b_new; + if(b_new > 0) + loss_new += C[GETI(ind)]*b_new*b_new; + x++; + } + } + + cond = cond + loss_new - loss_old; + if(cond <= 0) + break; + else + { + d_old = d; + d *= 0.5; + delta *= 0.5; + } + } + + w[j] += d; + + // recompute b[] if line search takes too many steps + if(num_linesearch >= max_num_linesearch) + { + info("#"); + for(int i=0; ix[i]; + while(x->index != -1) + { + b[x->index-1] -= w[i]*x->value; + x++; + } + } + } + } + + if(iter == 0) + Gnorm1_init = Gnorm1_new; + iter++; + if(iter % 10 == 0) + info("."); + + if(Gnorm1_new <= eps*Gnorm1_init) + { + if(active_size == w_size) + break; + else + { + active_size = w_size; + info("*"); + Gmax_old = INF; + continue; + } + } + + Gmax_old = Gmax_new; + } + + info("\noptimization finished, #iter = %d\n", iter); + if(iter >= max_iter) + info("\nWARNING: reaching max number of iterations\n"); + + // calculate objective value + + double v = 0; + int nnz = 0; + for(j=0; jx[j]; + while(x->index != -1) + { + x->value *= prob_col->y[x->index-1]; // restore x->value + x++; + } + if(w[j] != 0) + { + v += fabs(w[j]); + nnz++; + } + } + for(j=0; j 0) + v += C[GETI(j)]*b[j]*b[j]; + + info("Objective value = %lf\n", v); + info("#nonzeros/#features = %d/%d\n", nnz, w_size); + + delete [] index; + delete [] y; + delete [] b; + delete [] xj_sq; +} + +// A coordinate descent algorithm for +// L1-regularized logistic regression problems +// +// min_w \sum |wj| + C \sum log(1+exp(-yi w^T xi)), +// +// Given: +// x, y, Cp, Cn +// eps is the stopping tolerance +// +// solution will be put in w +// +// See Yuan et al. (2011) and appendix of LIBLINEAR paper, Fan et al. (2008) + +#undef GETI +#define GETI(i) (y[i]+1) +// To support weights for instances, use GETI(i) (i) + +static void solve_l1r_lr( + const problem *prob_col, double *w, double eps, + double Cp, double Cn) +{ + int l = prob_col->l; + int w_size = prob_col->n; + int j, s, newton_iter=0, iter=0; + int max_newton_iter = 100; + int max_iter = 1000; + int max_num_linesearch = 20; + int active_size; + int QP_active_size; + + double nu = 1e-12; + double inner_eps = 1; + double sigma = 0.01; + double w_norm, w_norm_new; + double z, G, H; + double Gnorm1_init; + double Gmax_old = INF; + double Gmax_new, Gnorm1_new; + double QP_Gmax_old = INF; + double QP_Gmax_new, QP_Gnorm1_new; + double delta, negsum_xTd, cond; + + int *index = new int[w_size]; + schar *y = new schar[l]; + double *Hdiag = new double[w_size]; + double *Grad = new double[w_size]; + double *wpd = new double[w_size]; + double *xjneg_sum = new double[w_size]; + double *xTd = new double[l]; + double *exp_wTx = new double[l]; + double *exp_wTx_new = new double[l]; + double *tau = new double[l]; + double *D = new double[l]; + feature_node *x; + + double C[3] = {Cn,0,Cp}; + + // Initial w can be set here. + for(j=0; jy[j] > 0) + y[j] = 1; + else + y[j] = -1; + + exp_wTx[j] = 0; + } + + w_norm = 0; + for(j=0; jx[j]; + while(x->index != -1) + { + int ind = x->index-1; + double val = x->value; + exp_wTx[ind] += w[j]*val; + if(y[ind] == -1) + xjneg_sum[j] += C[GETI(ind)]*val; + x++; + } + } + for(j=0; jx[j]; + while(x->index != -1) + { + int ind = x->index-1; + Hdiag[j] += x->value*x->value*D[ind]; + tmp += x->value*tau[ind]; + x++; + } + Grad[j] = -tmp + xjneg_sum[j]; + + double Gp = Grad[j]+1; + double Gn = Grad[j]-1; + double violation = 0; + if(w[j] == 0) + { + if(Gp < 0) + violation = -Gp; + else if(Gn > 0) + violation = Gn; + //outer-level shrinking + else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) + { + active_size--; + swap(index[s], index[active_size]); + s--; + continue; + } + } + else if(w[j] > 0) + violation = fabs(Gp); + else + violation = fabs(Gn); + + Gmax_new = max(Gmax_new, violation); + Gnorm1_new += violation; + } + + if(newton_iter == 0) + Gnorm1_init = Gnorm1_new; + + if(Gnorm1_new <= eps*Gnorm1_init) + break; + + iter = 0; + QP_Gmax_old = INF; + QP_active_size = active_size; + + for(int i=0; ix[j]; + G = Grad[j] + (wpd[j]-w[j])*nu; + while(x->index != -1) + { + int ind = x->index-1; + G += x->value*D[ind]*xTd[ind]; + x++; + } + + double Gp = G+1; + double Gn = G-1; + double violation = 0; + if(wpd[j] == 0) + { + if(Gp < 0) + violation = -Gp; + else if(Gn > 0) + violation = Gn; + //inner-level shrinking + else if(Gp>QP_Gmax_old/l && Gn<-QP_Gmax_old/l) + { + QP_active_size--; + swap(index[s], index[QP_active_size]); + s--; + continue; + } + } + else if(wpd[j] > 0) + violation = fabs(Gp); + else + violation = fabs(Gn); + + QP_Gmax_new = max(QP_Gmax_new, violation); + QP_Gnorm1_new += violation; + + // obtain solution of one-variable problem + if(Gp < H*wpd[j]) + z = -Gp/H; + else if(Gn > H*wpd[j]) + z = -Gn/H; + else + z = -wpd[j]; + + if(fabs(z) < 1.0e-12) + continue; + z = min(max(z,-10.0),10.0); + + wpd[j] += z; + + x = prob_col->x[j]; + while(x->index != -1) + { + int ind = x->index-1; + xTd[ind] += x->value*z; + x++; + } + } + + iter++; + + if(QP_Gnorm1_new <= inner_eps*Gnorm1_init) + { + //inner stopping + if(QP_active_size == active_size) + break; + //active set reactivation + else + { + QP_active_size = active_size; + QP_Gmax_old = INF; + continue; + } + } + + QP_Gmax_old = QP_Gmax_new; + } + + if(iter >= max_iter) + info("WARNING: reaching max number of inner iterations\n"); + + delta = 0; + w_norm_new = 0; + for(j=0; j= max_num_linesearch) + { + for(int i=0; ix[i]; + while(x->index != -1) + { + exp_wTx[x->index-1] += w[i]*x->value; + x++; + } + } + + for(int i=0; i= max_newton_iter) + info("WARNING: reaching max number of iterations\n"); + + // calculate objective value + + double v = 0; + int nnz = 0; + for(j=0; jl; + int n = prob->n; + int nnz = 0; + int *col_ptr = new int[n+1]; + feature_node *x_space; + prob_col->l = l; + prob_col->n = n; + prob_col->y = new double[l]; + prob_col->x = new feature_node*[n]; + + for(i=0; iy[i] = prob->y[i]; + + for(i=0; ix[i]; + while(x->index != -1) + { + nnz++; + col_ptr[x->index]++; + x++; + } + } + for(i=1; ix[i] = &x_space[col_ptr[i]]; + + for(i=0; ix[i]; + while(x->index != -1) + { + int ind = x->index-1; + x_space[col_ptr[ind]].index = i+1; // starts from 1 + x_space[col_ptr[ind]].value = x->value; + col_ptr[ind]++; + x++; + } + } + for(i=0; il; + int max_nr_class = 16; + int nr_class = 0; + int *label = Malloc(int,max_nr_class); + int *count = Malloc(int,max_nr_class); + int *data_label = Malloc(int,l); + int i; + + for(i=0;iy[i]; + int j; + for(j=0;jeps; + int pos = 0; + int neg = 0; + for(int i=0;il;i++) + if(prob->y[i] > 0) + pos++; + neg = prob->l - pos; + + double primal_solver_tol = eps*max(min(pos,neg), 1)/prob->l; + + function *fun_obj=NULL; + switch(param->solver_type) + { + case L2R_LR: + { + double *C = new double[prob->l]; + for(int i = 0; i < prob->l; i++) + { + if(prob->y[i] > 0) + C[i] = Cp; + else + C[i] = Cn; + } + fun_obj=new l2r_lr_fun(prob, C); + TRON tron_obj(fun_obj, primal_solver_tol); + tron_obj.set_print_string(liblinear_print_string); + tron_obj.tron(w); + delete fun_obj; + delete C; + break; + } + case L2R_L2LOSS_SVC: + { + double *C = new double[prob->l]; + for(int i = 0; i < prob->l; i++) + { + if(prob->y[i] > 0) + C[i] = Cp; + else + C[i] = Cn; + } + fun_obj=new l2r_l2_svc_fun(prob, C); + TRON tron_obj(fun_obj, primal_solver_tol); + tron_obj.set_print_string(liblinear_print_string); + tron_obj.tron(w); + delete fun_obj; + delete C; + break; + } + case L2R_L2LOSS_SVC_DUAL: + solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L2LOSS_SVC_DUAL); + break; + case L2R_L1LOSS_SVC_DUAL: + solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L1LOSS_SVC_DUAL); + break; + case L1R_L2LOSS_SVC: + { + problem prob_col; + feature_node *x_space = NULL; + transpose(prob, &x_space ,&prob_col); + solve_l1r_l2_svc(&prob_col, w, primal_solver_tol, Cp, Cn); + delete [] prob_col.y; + delete [] prob_col.x; + delete [] x_space; + break; + } + case L1R_LR: + { + problem prob_col; + feature_node *x_space = NULL; + transpose(prob, &x_space ,&prob_col); + solve_l1r_lr(&prob_col, w, primal_solver_tol, Cp, Cn); + delete [] prob_col.y; + delete [] prob_col.x; + delete [] x_space; + break; + } + case L2R_LR_DUAL: + solve_l2r_lr_dual(prob, w, eps, Cp, Cn); + break; + case L2R_L2LOSS_SVR: + { + double *C = new double[prob->l]; + for(int i = 0; i < prob->l; i++) + C[i] = param->C; + + fun_obj=new l2r_l2_svr_fun(prob, C, param->p); + TRON tron_obj(fun_obj, param->eps); + tron_obj.set_print_string(liblinear_print_string); + tron_obj.tron(w); + delete fun_obj; + delete C; + break; + + } + case L2R_L1LOSS_SVR_DUAL: + solve_l2r_l1l2_svr(prob, w, param, L2R_L1LOSS_SVR_DUAL); + break; + case L2R_L2LOSS_SVR_DUAL: + solve_l2r_l1l2_svr(prob, w, param, L2R_L2LOSS_SVR_DUAL); + break; + default: + fprintf(stderr, "ERROR: unknown solver_type\n"); + break; + } +} + +// +// Interface functions +// +model* train(const problem *prob, const parameter *param) +{ + int i,j; + int l = prob->l; + int n = prob->n; + int w_size = prob->n; + model *model_ = Malloc(model,1); + + if(prob->bias>=0) + model_->nr_feature=n-1; + else + model_->nr_feature=n; + model_->param = *param; + model_->bias = prob->bias; + + if(param->solver_type == L2R_L2LOSS_SVR || + param->solver_type == L2R_L1LOSS_SVR_DUAL || + param->solver_type == L2R_L2LOSS_SVR_DUAL) + { + model_->w = Malloc(double, w_size); + model_->nr_class = 2; + model_->label = NULL; + train_one(prob, param, &model_->w[0], 0, 0); + } + else + { + int nr_class; + int *label = NULL; + int *start = NULL; + int *count = NULL; + int *perm = Malloc(int,l); + + // group training data of the same class + group_classes(prob,&nr_class,&label,&start,&count,perm); + + model_->nr_class=nr_class; + model_->label = Malloc(int,nr_class); + for(i=0;ilabel[i] = label[i]; + + // calculate weighted C + double *weighted_C = Malloc(double, nr_class); + for(i=0;iC; + for(i=0;inr_weight;i++) + { + for(j=0;jweight_label[i] == label[j]) + break; + if(j == nr_class) + fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]); + else + weighted_C[j] *= param->weight[i]; + } + + // constructing the subproblem + feature_node **x = Malloc(feature_node *,l); + for(i=0;ix[perm[i]]; + + int k; + problem sub_prob; + sub_prob.l = l; + sub_prob.n = n; + sub_prob.x = Malloc(feature_node *,sub_prob.l); + sub_prob.y = Malloc(double,sub_prob.l); + + for(k=0; ksolver_type == MCSVM_CS) + { + model_->w=Malloc(double, n*nr_class); + for(i=0;ieps); + Solver.Solve(model_->w); + } + else + { + if(nr_class == 2) + { + model_->w=Malloc(double, w_size); + + int e0 = start[0]+count[0]; + k=0; + for(; kw[0], weighted_C[0], weighted_C[1]); + } + else + { + model_->w=Malloc(double, w_size*nr_class); + double *w=Malloc(double, w_size); + for(i=0;iC); + + for(int j=0;jw[j*nr_class+i] = w[j]; + } + free(w); + } + + } + + free(x); + free(label); + free(start); + free(count); + free(perm); + free(sub_prob.x); + free(sub_prob.y); + free(weighted_C); + } + return model_; +} + +void cross_validation(const problem *prob, const parameter *param, int nr_fold, double *target) +{ + int i; + int *fold_start = Malloc(int,nr_fold+1); + int l = prob->l; + int *perm = Malloc(int,l); + + for(i=0;ibias; + subprob.n = prob->n; + subprob.l = l-(end-begin); + subprob.x = Malloc(struct feature_node*,subprob.l); + subprob.y = Malloc(double,subprob.l); + + k=0; + for(j=0;jx[perm[j]]; + subprob.y[k] = prob->y[perm[j]]; + ++k; + } + for(j=end;jx[perm[j]]; + subprob.y[k] = prob->y[perm[j]]; + ++k; + } + struct model *submodel = train(&subprob,param); + for(j=begin;jx[perm[j]]); + free_and_destroy_model(&submodel); + free(subprob.x); + free(subprob.y); + } + free(fold_start); + free(perm); +} + +double predict_values(const struct model *model_, const struct feature_node *x, double *dec_values) +{ + int idx; + int n; + if(model_->bias>=0) + n=model_->nr_feature+1; + else + n=model_->nr_feature; + double *w=model_->w; + int nr_class=model_->nr_class; + int i; + int nr_w; + if(nr_class==2 && model_->param.solver_type != MCSVM_CS) + nr_w = 1; + else + nr_w = nr_class; + + const feature_node *lx=x; + for(i=0;iindex)!=-1; lx++) + { + // the dimension of testing data may exceed that of training + if(idx<=n) + for(i=0;ivalue; + } + + if(nr_class==2) + { + if(model_->param.solver_type == L2R_L2LOSS_SVR || + model_->param.solver_type == L2R_L1LOSS_SVR_DUAL || + model_->param.solver_type == L2R_L2LOSS_SVR_DUAL) + return dec_values[0]; + else + return (dec_values[0]>0)?model_->label[0]:model_->label[1]; + } + else + { + int dec_max_idx = 0; + for(i=1;i dec_values[dec_max_idx]) + dec_max_idx = i; + } + return model_->label[dec_max_idx]; + } +} + +double predict(const model *model_, const feature_node *x) +{ + double *dec_values = Malloc(double, model_->nr_class); + double label=predict_values(model_, x, dec_values); + free(dec_values); + return label; +} + +double predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates) +{ + if(check_probability_model(model_)) + { + int i; + int nr_class=model_->nr_class; + int nr_w; + if(nr_class==2) + nr_w = 1; + else + nr_w = nr_class; + + double label=predict_values(model_, x, prob_estimates); + for(i=0;inr_feature; + int n; + const parameter& param = model_->param; + + if(model_->bias>=0) + n=nr_feature+1; + else + n=nr_feature; + int w_size = n; + FILE *fp = fopen(model_file_name,"w"); + if(fp==NULL) return -1; + + char *old_locale = strdup(setlocale(LC_ALL, NULL)); + setlocale(LC_ALL, "C"); + + int nr_w; + if(model_->nr_class==2 && model_->param.solver_type != MCSVM_CS) + nr_w=1; + else + nr_w=model_->nr_class; + + fprintf(fp, "solver_type %s\n", solver_type_table[param.solver_type]); + fprintf(fp, "nr_class %d\n", model_->nr_class); + + if(model_->label) + { + fprintf(fp, "label"); + for(i=0; inr_class; i++) + fprintf(fp, " %d", model_->label[i]); + fprintf(fp, "\n"); + } + + fprintf(fp, "nr_feature %d\n", nr_feature); + + fprintf(fp, "bias %.16g\n", model_->bias); + + fprintf(fp, "w\n"); + for(i=0; iw[i*nr_w+j]); + fprintf(fp, "\n"); + } + + setlocale(LC_ALL, old_locale); + free(old_locale); + + if (ferror(fp) != 0 || fclose(fp) != 0) return -1; + else return 0; +} + +struct model *load_model(const char *model_file_name) +{ + FILE *fp = fopen(model_file_name,"r"); + if(fp==NULL) return NULL; + + int i; + int nr_feature; + int n; + int nr_class; + double bias; + model *model_ = Malloc(model,1); + parameter& param = model_->param; + + model_->label = NULL; + + char *old_locale = strdup(setlocale(LC_ALL, NULL)); + setlocale(LC_ALL, "C"); + + char cmd[81]; + while(1) + { + fscanf(fp,"%80s",cmd); + if(strcmp(cmd,"solver_type")==0) + { + fscanf(fp,"%80s",cmd); + int i; + for(i=0;solver_type_table[i];i++) + { + if(strcmp(solver_type_table[i],cmd)==0) + { + param.solver_type=i; + break; + } + } + if(solver_type_table[i] == NULL) + { + fprintf(stderr,"unknown solver type.\n"); + + setlocale(LC_ALL, old_locale); + free(model_->label); + free(model_); + free(old_locale); + return NULL; + } + } + else if(strcmp(cmd,"nr_class")==0) + { + fscanf(fp,"%d",&nr_class); + model_->nr_class=nr_class; + } + else if(strcmp(cmd,"nr_feature")==0) + { + fscanf(fp,"%d",&nr_feature); + model_->nr_feature=nr_feature; + } + else if(strcmp(cmd,"bias")==0) + { + fscanf(fp,"%lf",&bias); + model_->bias=bias; + } + else if(strcmp(cmd,"w")==0) + { + break; + } + else if(strcmp(cmd,"label")==0) + { + int nr_class = model_->nr_class; + model_->label = Malloc(int,nr_class); + for(int i=0;ilabel[i]); + } + else + { + fprintf(stderr,"unknown text in model file: [%s]\n",cmd); + setlocale(LC_ALL, old_locale); + free(model_->label); + free(model_); + free(old_locale); + return NULL; + } + } + + nr_feature=model_->nr_feature; + if(model_->bias>=0) + n=nr_feature+1; + else + n=nr_feature; + int w_size = n; + int nr_w; + if(nr_class==2 && param.solver_type != MCSVM_CS) + nr_w = 1; + else + nr_w = nr_class; + + model_->w=Malloc(double, w_size*nr_w); + for(i=0; iw[i*nr_w+j]); + fscanf(fp, "\n"); + } + + setlocale(LC_ALL, old_locale); + free(old_locale); + + if (ferror(fp) != 0 || fclose(fp) != 0) return NULL; + + return model_; +} + +int get_nr_feature(const model *model_) +{ + return model_->nr_feature; +} + +int get_nr_class(const model *model_) +{ + return model_->nr_class; +} + +void get_labels(const model *model_, int* label) +{ + if (model_->label != NULL) + for(int i=0;inr_class;i++) + label[i] = model_->label[i]; +} + +void free_model_content(struct model *model_ptr) +{ + if(model_ptr->w != NULL) + free(model_ptr->w); + if(model_ptr->label != NULL) + free(model_ptr->label); +} + +void free_and_destroy_model(struct model **model_ptr_ptr) +{ + struct model *model_ptr = *model_ptr_ptr; + if(model_ptr != NULL) + { + free_model_content(model_ptr); + free(model_ptr); + } +} + +void destroy_param(parameter* param) +{ + if(param->weight_label != NULL) + free(param->weight_label); + if(param->weight != NULL) + free(param->weight); +} + +const char *check_parameter(const problem *prob, const parameter *param) +{ + if(param->eps <= 0) + return "eps <= 0"; + + if(param->C <= 0) + return "C <= 0"; + + if(param->p < 0) + return "p < 0"; + + if(param->solver_type != L2R_LR + && param->solver_type != L2R_L2LOSS_SVC_DUAL + && param->solver_type != L2R_L2LOSS_SVC + && param->solver_type != L2R_L1LOSS_SVC_DUAL + && param->solver_type != MCSVM_CS + && param->solver_type != L1R_L2LOSS_SVC + && param->solver_type != L1R_LR + && param->solver_type != L2R_LR_DUAL + && param->solver_type != L2R_L2LOSS_SVR + && param->solver_type != L2R_L2LOSS_SVR_DUAL + && param->solver_type != L2R_L1LOSS_SVR_DUAL) + return "unknown solver type"; + + return NULL; +} + +int check_probability_model(const struct model *model_) +{ + return (model_->param.solver_type==L2R_LR || + model_->param.solver_type==L2R_LR_DUAL || + model_->param.solver_type==L1R_LR); +} + +void set_print_string_function(void (*print_func)(const char*)) +{ + if (print_func == NULL) + liblinear_print_string = &print_string_stdout; + else + liblinear_print_string = print_func; +} + diff --git a/random/include/LibLinear/linear.h b/random/include/LibLinear/linear.h new file mode 100755 index 0000000..f02a93e --- /dev/null +++ b/random/include/LibLinear/linear.h @@ -0,0 +1,74 @@ +#ifndef _LIBLINEAR_H +#define _LIBLINEAR_H + +#ifdef __cplusplus +extern "C" { +#endif + +struct feature_node +{ + int index; + double value; +}; + +struct problem +{ + int l, n; + double *y; + struct feature_node **x; + double bias; /* < 0 if no bias term */ +}; + +enum { L2R_LR, L2R_L2LOSS_SVC_DUAL, L2R_L2LOSS_SVC, L2R_L1LOSS_SVC_DUAL, MCSVM_CS, L1R_L2LOSS_SVC, L1R_LR, L2R_LR_DUAL, L2R_L2LOSS_SVR = 11, L2R_L2LOSS_SVR_DUAL, L2R_L1LOSS_SVR_DUAL }; /* solver_type */ + +struct parameter +{ + int solver_type; + + /* these are for training only */ + double eps; /* stopping criteria */ + double C; + int nr_weight; + int *weight_label; + double* weight; + double p; +}; + +struct model +{ + struct parameter param; + int nr_class; /* number of classes */ + int nr_feature; + double *w; + int *label; /* label of each class */ + double bias; +}; + +struct model* train(const struct problem *prob, const struct parameter *param); +void cross_validation(const struct problem *prob, const struct parameter *param, int nr_fold, double *target); + +double predict_values(const struct model *model_, const struct feature_node *x, double* dec_values); +double predict(const struct model *model_, const struct feature_node *x); +double predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates); + +int save_model(const char *model_file_name, const struct model *model_); +struct model *load_model(const char *model_file_name); + +int get_nr_feature(const struct model *model_); +int get_nr_class(const struct model *model_); +void get_labels(const struct model *model_, int* label); + +void free_model_content(struct model *model_ptr); +void free_and_destroy_model(struct model **model_ptr_ptr); +void destroy_param(struct parameter *param); + +const char *check_parameter(const struct problem *prob, const struct parameter *param); +int check_probability_model(const struct model *model); +void set_print_string_function(void (*print_func) (const char*)); + +#ifdef __cplusplus +} +#endif + +#endif /* _LIBLINEAR_H */ + diff --git a/random/include/LibLinear/train.c b/random/include/LibLinear/train.c new file mode 100755 index 0000000..27bcc9c --- /dev/null +++ b/random/include/LibLinear/train.c @@ -0,0 +1,401 @@ +#include +#include +#include +#include +#include +#include +#include "linear.h" +#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) +#define INF HUGE_VAL + +#pragma warning(disable:4996) + +void print_null(const char *s) {} + +void exit_with_help() +{ + printf( + "Usage: train [options] training_set_file [model_file]\n" + "options:\n" + "-s type : set type of solver (default 1)\n" + " for multi-class classification\n" + " 0 -- L2-regularized logistic regression (primal)\n" + " 1 -- L2-regularized L2-loss support vector classification (dual)\n" + " 2 -- L2-regularized L2-loss support vector classification (primal)\n" + " 3 -- L2-regularized L1-loss support vector classification (dual)\n" + " 4 -- support vector classification by Crammer and Singer\n" + " 5 -- L1-regularized L2-loss support vector classification\n" + " 6 -- L1-regularized logistic regression\n" + " 7 -- L2-regularized logistic regression (dual)\n" + " for regression\n" + " 11 -- L2-regularized L2-loss support vector regression (primal)\n" + " 12 -- L2-regularized L2-loss support vector regression (dual)\n" + " 13 -- L2-regularized L1-loss support vector regression (dual)\n" + "-c cost : set the parameter C (default 1)\n" + "-p epsilon : set the epsilon in loss function of SVR (default 0.1)\n" + "-e epsilon : set tolerance of termination criterion\n" + " -s 0 and 2\n" + " |f'(w)|_2 <= eps*min(pos,neg)/l*|f'(w0)|_2,\n" + " where f is the primal function and pos/neg are # of\n" + " positive/negative data (default 0.01)\n" + " -s 11\n" + " |f'(w)|_2 <= eps*|f'(w0)|_2 (default 0.001)\n" + " -s 1, 3, 4, and 7\n" + " Dual maximal violation <= eps; similar to libsvm (default 0.1)\n" + " -s 5 and 6\n" + " |f'(w)|_1 <= eps*min(pos,neg)/l*|f'(w0)|_1,\n" + " where f is the primal function (default 0.01)\n" + " -s 12 and 13\n" + " |f'(alpha)|_1 <= eps |f'(alpha0)|,\n" + " where f is the dual function (default 0.1)\n" + "-B bias : if bias >= 0, instance x becomes [x; bias]; if < 0, no bias term added (default -1)\n" + "-wi weight: weights adjust the parameter C of different classes (see README for details)\n" + "-v n: n-fold cross validation mode\n" + "-q : quiet mode (no outputs)\n" + ); + exit(1); +} + +void exit_input_error(int line_num) +{ + fprintf(stderr,"Wrong input format at line %d\n", line_num); + exit(1); +} + +static char *line = NULL; +static int max_line_len; + +static char* readline(FILE *input) +{ + int len; + + if(fgets(line,max_line_len,input) == NULL) + return NULL; + + while(strrchr(line,'\n') == NULL) + { + max_line_len *= 2; + line = (char *) realloc(line,max_line_len); + len = (int) strlen(line); + if(fgets(line+len,max_line_len-len,input) == NULL) + break; + } + return line; +} + +void parse_command_line(int argc, char **argv, char *input_file_name, char *model_file_name); +void read_problem(const char *filename); +void do_cross_validation(); + +struct feature_node *x_space; +struct parameter param; +struct problem prob; +struct model* model_; +int flag_cross_validation; +int nr_fold; +double bias; + +int main(int argc, char **argv) +{ + char input_file_name[1024]; + char model_file_name[1024]; + const char *error_msg; + + parse_command_line(argc, argv, input_file_name, model_file_name); + read_problem(input_file_name); + error_msg = check_parameter(&prob,¶m); + + if(error_msg) + { + fprintf(stderr,"ERROR: %s\n",error_msg); + exit(1); + } + + if(flag_cross_validation) + { + do_cross_validation(); + } + else + { + model_=train(&prob, ¶m); + if(save_model(model_file_name, model_)) + { + fprintf(stderr,"can't save model to file %s\n",model_file_name); + exit(1); + } + free_and_destroy_model(&model_); + } + destroy_param(¶m); + free(prob.y); + free(prob.x); + free(x_space); + free(line); + + return 0; +} + +void do_cross_validation() +{ + int i; + int total_correct = 0; + double total_error = 0; + double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0; + double *target = Malloc(double, prob.l); + + cross_validation(&prob,¶m,nr_fold,target); + if(param.solver_type == L2R_L2LOSS_SVR || + param.solver_type == L2R_L1LOSS_SVR_DUAL || + param.solver_type == L2R_L2LOSS_SVR_DUAL) + { + for(i=0;i=argc) + exit_with_help(); + switch(argv[i-1][1]) + { + case 's': + param.solver_type = atoi(argv[i]); + break; + + case 'c': + param.C = atof(argv[i]); + break; + + case 'p': + param.p = atof(argv[i]); + break; + + case 'e': + param.eps = atof(argv[i]); + break; + + case 'B': + bias = atof(argv[i]); + break; + + case 'w': + ++param.nr_weight; + param.weight_label = (int *) realloc(param.weight_label,sizeof(int)*param.nr_weight); + param.weight = (double *) realloc(param.weight,sizeof(double)*param.nr_weight); + param.weight_label[param.nr_weight-1] = atoi(&argv[i-1][2]); + param.weight[param.nr_weight-1] = atof(argv[i]); + break; + + case 'v': + flag_cross_validation = 1; + nr_fold = atoi(argv[i]); + if(nr_fold < 2) + { + fprintf(stderr,"n-fold cross validation: n must >= 2\n"); + exit_with_help(); + } + break; + + case 'q': + print_func = &print_null; + i--; + break; + + default: + fprintf(stderr,"unknown option: -%c\n", argv[i-1][1]); + exit_with_help(); + break; + } + } + + set_print_string_function(print_func); + + // determine filenames + if(i>=argc) + exit_with_help(); + + strcpy(input_file_name, argv[i]); + + if(i max_index) + max_index = inst_max_index; + + if(prob.bias >= 0) + x_space[j++].value = prob.bias; + + x_space[j++].index = -1; + } + + if(prob.bias >= 0) + { + prob.n=max_index+1; + for(i=1;iindex = prob.n; + x_space[j-2].index = prob.n; + } + else + prob.n=max_index; + + fclose(fp); +} diff --git a/random/include/LibLinear/tron.cpp b/random/include/LibLinear/tron.cpp new file mode 100755 index 0000000..c678990 --- /dev/null +++ b/random/include/LibLinear/tron.cpp @@ -0,0 +1,235 @@ +#include +#include +#include +#include +#include "tron.h" + +#ifndef min +template static inline T min(T x,T y) { return (x static inline T max(T x,T y) { return (x>y)?x:y; } +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +extern double dnrm2_(int *, double *, int *); +extern double ddot_(int *, double *, int *, double *, int *); +extern int daxpy_(int *, double *, double *, int *, double *, int *); +extern int dscal_(int *, double *, double *, int *); + +#ifdef __cplusplus +} +#endif + +static void default_print(const char *buf) +{ + fputs(buf,stdout); + fflush(stdout); +} + +void TRON::info(const char *fmt,...) +{ + char buf[BUFSIZ]; + va_list ap; + va_start(ap,fmt); + vsprintf(buf,fmt,ap); + va_end(ap); + (*tron_print_string)(buf); +} + +TRON::TRON(const function *fun_obj, double eps, int max_iter) +{ + this->fun_obj=const_cast(fun_obj); + this->eps=eps; + this->max_iter=max_iter; + tron_print_string = default_print; +} + +TRON::~TRON() +{ +} + +void TRON::tron(double *w) +{ + // Parameters for updating the iterates. + double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75; + + // Parameters for updating the trust region size delta. + double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4; + + int n = fun_obj->get_nr_variable(); + int i, cg_iter; + double delta, snorm, one=1.0; + double alpha, f, fnew, prered, actred, gs; + int search = 1, iter = 1, inc = 1; + double *s = new double[n]; + double *r = new double[n]; + double *w_new = new double[n]; + double *g = new double[n]; + + for (i=0; ifun(w); + fun_obj->grad(w, g); + delta = dnrm2_(&n, g, &inc); + double gnorm1 = delta; + double gnorm = gnorm1; + + if (gnorm <= eps*gnorm1) + search = 0; + + iter = 1; + + while (iter <= max_iter && search) + { + cg_iter = trcg(delta, g, s, r); + + memcpy(w_new, w, sizeof(double)*n); + daxpy_(&n, &one, s, &inc, w_new, &inc); + + gs = ddot_(&n, g, &inc, s, &inc); + prered = -0.5*(gs-ddot_(&n, s, &inc, r, &inc)); + fnew = fun_obj->fun(w_new); + + // Compute the actual reduction. + actred = f - fnew; + + // On the first iteration, adjust the initial step bound. + snorm = dnrm2_(&n, s, &inc); + if (iter == 1) + delta = min(delta, snorm); + + // Compute prediction alpha*snorm of the step. + if (fnew - f - gs <= 0) + alpha = sigma3; + else + alpha = max(sigma1, -0.5*(gs/(fnew - f - gs))); + + // Update the trust region bound according to the ratio of actual to predicted reduction. + if (actred < eta0*prered) + delta = min(max(alpha, sigma1)*snorm, sigma2*delta); + else if (actred < eta1*prered) + delta = max(sigma1*delta, min(alpha*snorm, sigma2*delta)); + else if (actred < eta2*prered) + delta = max(sigma1*delta, min(alpha*snorm, sigma3*delta)); + else + delta = max(delta, min(alpha*snorm, sigma3*delta)); + + info("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter); + + if (actred > eta0*prered) + { + iter++; + memcpy(w, w_new, sizeof(double)*n); + f = fnew; + fun_obj->grad(w, g); + + gnorm = dnrm2_(&n, g, &inc); + if (gnorm <= eps*gnorm1) + break; + } + if (f < -1.0e+32) + { + info("WARNING: f < -1.0e+32\n"); + break; + } + if (fabs(actred) <= 0 && prered <= 0) + { + info("WARNING: actred and prered <= 0\n"); + break; + } + if (fabs(actred) <= 1.0e-12*fabs(f) && + fabs(prered) <= 1.0e-12*fabs(f)) + { + info("WARNING: actred and prered too small\n"); + break; + } + } + + delete[] g; + delete[] r; + delete[] w_new; + delete[] s; +} + +int TRON::trcg(double delta, double *g, double *s, double *r) +{ + int i, inc = 1; + int n = fun_obj->get_nr_variable(); + double one = 1; + double *d = new double[n]; + double *Hd = new double[n]; + double rTr, rnewTrnew, alpha, beta, cgtol; + + for (i=0; iHv(d, Hd); + + alpha = rTr/ddot_(&n, d, &inc, Hd, &inc); + daxpy_(&n, &alpha, d, &inc, s, &inc); + if (dnrm2_(&n, s, &inc) > delta) + { + info("cg reaches trust region boundary\n"); + alpha = -alpha; + daxpy_(&n, &alpha, d, &inc, s, &inc); + + double std = ddot_(&n, s, &inc, d, &inc); + double sts = ddot_(&n, s, &inc, s, &inc); + double dtd = ddot_(&n, d, &inc, d, &inc); + double dsq = delta*delta; + double rad = sqrt(std*std + dtd*(dsq-sts)); + if (std >= 0) + alpha = (dsq - sts)/(std + rad); + else + alpha = (rad - std)/dtd; + daxpy_(&n, &alpha, d, &inc, s, &inc); + alpha = -alpha; + daxpy_(&n, &alpha, Hd, &inc, r, &inc); + break; + } + alpha = -alpha; + daxpy_(&n, &alpha, Hd, &inc, r, &inc); + rnewTrnew = ddot_(&n, r, &inc, r, &inc); + beta = rnewTrnew/rTr; + dscal_(&n, &beta, d, &inc); + daxpy_(&n, &one, r, &inc, d, &inc); + rTr = rnewTrnew; + } + + delete[] d; + delete[] Hd; + + return(cg_iter); +} + +double TRON::norm_inf(int n, double *x) +{ + double dmax = fabs(x[0]); + for (int i=1; i= dmax) + dmax = fabs(x[i]); + return(dmax); +} + +void TRON::set_print_string(void (*print_string) (const char *buf)) +{ + tron_print_string = print_string; +} diff --git a/random/include/LibLinear/tron.h b/random/include/LibLinear/tron.h new file mode 100755 index 0000000..879e0f1 --- /dev/null +++ b/random/include/LibLinear/tron.h @@ -0,0 +1,36 @@ +#ifndef _TRON_H +#define _TRON_H + +#pragma warning(disable:4996) + +class function +{ +public: + virtual double fun(double *w) = 0 ; + virtual void grad(double *w, double *g) = 0 ; + virtual void Hv(double *s, double *Hs) = 0 ; + + virtual int get_nr_variable(void) = 0 ; + virtual ~function(void){} +}; + +class TRON +{ +public: + TRON(const function *fun_obj, double eps = 0.1, int max_iter = 1000); + ~TRON(); + + void tron(double *w); + void set_print_string(void (*i_print) (const char *buf)); + +private: + int trcg(double delta, double *g, double *s, double *r); + double norm_inf(int n, double *x); + + double eps; + int max_iter; + function *fun_obj; + void info(const char *fmt,...); + void (*tron_print_string)(const char *buf); +}; +#endif diff --git a/random/include/Objectness/CmFile.h b/random/include/Objectness/CmFile.h new file mode 100755 index 0000000..80c2578 --- /dev/null +++ b/random/include/Objectness/CmFile.h @@ -0,0 +1,170 @@ +#pragma once + +struct CmFile +{ + //static string BrowseFile(const char* strFilter = "Images (*.jpg;*.png)\0*.jpg;*.png\0All (*.*)\0*.*\0\0", bool isOpen = true); + //static string BrowseFolder(); + + static inline string GetFolder(CStr& path); + static inline string GetName(CStr& path); + static inline string GetNameNE(CStr& path); + static inline string GetPathNE(CStr& path); + + // Get file names from a wildcard. Eg: GetNames("D:\\*.jpg", imgNames); +#ifdef _WIN32 + static int GetNames(CStr &nameW, vecS &names, string &dir = string()); +#else + static int GetNames(CStr &nameW, vecS &names, string &dir); + static int GetNames(CStr &nameW, vecS &names); +#endif + static int GetNames(CStr& rootFolder, CStr &fileW, vecS &names); +#ifdef _WIN32 + static int GetNamesNE(CStr& nameWC, vecS &names, string &dir = string(), string &ext = string()); +#else + static int GetNamesNE(CStr& nameWC, vecS &names, string &dir, string &ext); + static int GetNamesNE(CStr& nameWC, vecS &names); +#endif + static int GetNamesNE(CStr& rootFolder, CStr &fileW, vecS &names); + static inline string GetExtention(CStr name); + + //static inline bool FileExist(CStr& filePath); + //static inline bool FilesExist(CStr& fileW); + //static inline bool FolderExist(CStr& strPath); + + static inline string GetWkDir(); + + static BOOL MkDir(CStr& path); + + // Eg: RenameImages("D:/DogImages/*.jpg", "F:/Images", "dog", ".jpg"); + //static int Rename(CStr& srcNames, CStr& dstDir, const char* nameCommon, const char* nameExt); + + //static inline void RmFile(CStr& fileW); + //static void RmFolder(CStr& dir); + //static void CleanFolder(CStr& dir, bool subFolder = false); + + static int GetSubFolders(CStr& folder, vecS& subFolders); + + //inline static bool Copy(CStr &src, CStr &dst, bool failIfExist = false); + //inline static BOOL Move(CStr &src, CStr &dst, DWORD dwFlags = MOVEFILE_REPLACE_EXISTING | MOVEFILE_COPY_ALLOWED | MOVEFILE_WRITE_THROUGH); + //static BOOL Move2Dir(CStr &srcW, CStr dstDir); + //static BOOL Copy2Dir(CStr &srcW, CStr dstDir); + + //Load mask image and threshold thus noisy by compression can be removed + //static Mat LoadMask(CStr& fileName); + + //static void WriteNullFile(CStr& fileName) {FILE *f = fopen(_S(fileName), "w"); fclose(f);} + + //static void ChkImgs(CStr &imgW); + + //static void RunProgram(CStr &fileName, CStr ¶meters = "", bool waiteF = false, bool showW = true); + + //static void SegOmpThrdNum(double ratio = 0.8); + + // Copy files and add suffix. e.g. copyAddSuffix("./*.jpg", "./Imgs/", "_Img.jpg") + //static void copyAddSuffix(CStr &srcW, CStr &dstDir, CStr &dstSuffix); + + static vecS loadStrList(CStr &fName); + static bool writeStrList(CStr &fName, const vecS &strs); +}; + +/************************************************************************/ +/* Implementation of inline functions */ +/************************************************************************/ +string CmFile::GetFolder(CStr& path) +{ + return path.substr(0, path.find_last_of("\\/")+1); +} + +string CmFile::GetName(CStr& path) +{ + int start = path.find_last_of("\\/")+1; + int end = path.find_last_not_of(' ')+1; + return path.substr(start, end - start); +} + +string CmFile::GetNameNE(CStr& path) +{ + int start = path.find_last_of("\\/")+1; + int end = path.find_last_of('.'); + if (end >= 0) + return path.substr(start, end - start); + else + return path.substr(start, path.find_last_not_of(' ')+1 - start); +} + +string CmFile::GetPathNE(CStr& path) +{ + int end = path.find_last_of('.'); + if (end >= 0) + return path.substr(0, end); + else + return path.substr(0, path.find_last_not_of(' ') + 1); +} + +string CmFile::GetExtention(CStr name) +{ + return name.substr(name.find_last_of('.')); +} + +/*BOOL CmFile::Copy(CStr &src, CStr &dst, bool failIfExist) +{ + return ::CopyFileA(src.c_str(), dst.c_str(), failIfExist); +}*/ + +/*BOOL CmFile::Move(CStr &src, CStr &dst, DWORD dwFlags) +{ + return MoveFileExA(src.c_str(), dst.c_str(), dwFlags); +}*/ + +/*void CmFile::RmFile(CStr& fileW) +{ + vecS names; + string dir; + int fNum = CmFile::GetNames(fileW, names, dir); + for (int i = 0; i < fNum; i++) + ::DeleteFileA(_S(dir + names[i])); +}*/ + + +// Test whether a file exist +/*bool CmFile::FileExist(CStr& filePath) +{ + if (filePath.size() == 0) + return false; + + return GetFileAttributesA(_S(filePath)) != INVALID_FILE_ATTRIBUTES; // || GetLastError() != ERROR_FILE_NOT_FOUND; +}*/ + +/*bool CmFile::FilesExist(CStr& fileW) +{ + vecS names; + int fNum = GetNames(fileW, names); + return fNum > 0; +}*/ + +/*string CmFile::GetWkDir() +{ + string wd; + wd.resize(1024); + DWORD len = GetCurrentDirectoryA(1024, &wd[0]); + wd.resize(len); + return wd; +}*/ + +/*bool CmFile::FolderExist(CStr& strPath) +{ + int i = (int)strPath.size() - 1; + for (; i >= 0 && (strPath[i] == '\\' || strPath[i] == '/'); i--) + ; + string str = strPath.substr(0, i+1); + + WIN32_FIND_DATAA wfd; + HANDLE hFind = FindFirstFileA(_S(str), &wfd); + bool rValue = (hFind != INVALID_HANDLE_VALUE) && (wfd.dwFileAttributes & FILE_ATTRIBUTE_DIRECTORY); + FindClose(hFind); + return rValue; +}*/ + +/************************************************************************/ +/* Implementations */ +/************************************************************************/ diff --git a/random/include/Objectness/CmShow.h b/random/include/Objectness/CmShow.h new file mode 100755 index 0000000..11d99a0 --- /dev/null +++ b/random/include/Objectness/CmShow.h @@ -0,0 +1,9 @@ +#pragma once +class CmShow +{ +public: + static Mat HistBins(CMat& color3f, CMat& val, CStr& title, bool descendShow = false, CMat &with = Mat()); + static void showTinyMat(CStr &title, CMat &m); + static inline void SaveShow(CMat& img, CStr& title); +}; + diff --git a/random/include/Objectness/CmTimer.h b/random/include/Objectness/CmTimer.h new file mode 100755 index 0000000..0fc9b41 --- /dev/null +++ b/random/include/Objectness/CmTimer.h @@ -0,0 +1,83 @@ +#pragma once + +class CmTimer +{ +public: + CmTimer(CStr t):title(t) { is_started = false; start_clock = 0; cumulative_clock = 0; n_starts = 0; } + + ~CmTimer(){ if (is_started) printf("CmTimer '%s' is started and is being destroyed.\n", title.c_str()); } + + inline void Start(); + inline void Stop(); + inline void Reset(); + + inline bool Report(); + inline bool StopAndReport() { Stop(); return Report(); } + inline float TimeInSeconds(); + +private: + CStr title; + + bool is_started; + clock_t start_clock; + clock_t cumulative_clock; + unsigned int n_starts; +}; + +/************************************************************************/ +/* Implementations */ +/************************************************************************/ + +void CmTimer::Start() +{ + if (is_started){ + printf("CmTimer '%s' is already started. Nothing done.\n", title.c_str()); + return; + } + + is_started = true; + n_starts++; + start_clock = clock(); +} + +void CmTimer::Stop() +{ + if (!is_started){ + printf("CmTimer '%s' is started. Nothing done\n", title.c_str()); + return; + } + + cumulative_clock += clock() - start_clock; + is_started = false; +} + +void CmTimer::Reset() +{ + if (is_started) { + printf("CmTimer '%s'is started during reset request.\n Only reset cumulative time.\n", title.c_str()); + return; + } + cumulative_clock = 0; +} + +bool CmTimer::Report() +{ + if (is_started){ + printf("CmTimer '%s' is started.\n Cannot provide a time report.", title.c_str()); + return false; + } + + float timeUsed = TimeInSeconds(); + printf("[%s] CumuTime: %gs, #run: %d, AvgTime: %gs\n", title.c_str(), timeUsed, n_starts, timeUsed/n_starts); + return true; +} + +float CmTimer::TimeInSeconds() +{ + if (is_started){ + printf("CmTimer '%s' is started. Nothing done\n", title.c_str()); + return 0; + } + return float(cumulative_clock) / CLOCKS_PER_SEC; +} + diff --git a/random/include/Objectness/DataSetVOC.h b/random/include/Objectness/DataSetVOC.h new file mode 100755 index 0000000..1f33539 --- /dev/null +++ b/random/include/Objectness/DataSetVOC.h @@ -0,0 +1,77 @@ +#pragma once + +struct DataSetVOC +{ + DataSetVOC(CStr &wkDir); + ~DataSetVOC(void); + + // Organization structure data for the dataset + string wkDir; // Root working directory, all other directories are relative to this one + string resDir, localDir; // Directory for saving results and local data + string imgPathW, annoPathW; // Image and annotation path + + // Information for training and testing + int trainNum, testNum; + vecS trainSet, testSet; // File names (NE) for training and testing images + vecS classNames; // Object class names + vector > gtTrainBoxes, gtTestBoxes; // Ground truth bounding boxes for training and testing images + vector gtTrainClsIdx, gtTestClsIdx; // Object class indexes + + + // Load annotations + void loadAnnotations(); + + static bool cvt2OpenCVYml(CStr &annoDir); // Needs to call yml.m in this solution before running this function. + + //static bool importSaliencyBench(CStr &salDir = "./THUS10000/", CStr &vocDir = "./THUS10000/"); + + //static void importPaulData(CStr &inData = "Z:/datasets/toMing/", CStr &outData = "D:/WkDir/DetectionProposals/Paul/"); + + //static bool importImageNetBenchMark(CStr &orgDir = "D:/WkDir/ImageNet/", CStr &newDir = "D:/WkDir/DetectionProposals/ImgNet/"); + + static inline double interUnio(const Vec4i &box1, const Vec4i &box2); + + // Get training and testing for demonstrating the generative of the objectness over classes + void getTrainTest(); + +public: // Used for testing the ability of generic over classes + void loadDataGenericOverCls(); + +private: + void loadBox(const FileNode &fn, vector &boxes, vecI &clsIdx); + bool loadBBoxes(CStr &nameNE, vector &boxes, vecI &clsIdx); + static void getXmlStrVOC(CStr &fName, string &buf); + static inline string keepXmlChar(CStr &str); + static bool cvt2OpenCVYml(CStr &yamlName, CStr &ymlName); // Needs to call yml.m in this solution before running this function. +}; + +string DataSetVOC::keepXmlChar(CStr &_str) +{ + string str = _str; + int sz = (int)str.size(), count = 0; + for (int i = 0; i < sz; i++){ + char c = str[i]; + if ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') || (c >= '0' && c <= '9') || c == ' ' || c == '.') + str[count++] = str[i]; + } + str.resize(count); + return str; +} + +double DataSetVOC::interUnio(const Vec4i &bb, const Vec4i &bbgt) +{ + int bi[4]; + bi[0] = max(bb[0], bbgt[0]); + bi[1] = max(bb[1], bbgt[1]); + bi[2] = min(bb[2], bbgt[2]); + bi[3] = min(bb[3], bbgt[3]); + + double iw = bi[2] - bi[0] + 1; + double ih = bi[3] - bi[1] + 1; + double ov = 0; + if (iw>0 && ih>0){ + double ua = (bb[2]-bb[0]+1)*(bb[3]-bb[1]+1)+(bbgt[2]-bbgt[0]+1)*(bbgt[3]-bbgt[1]+1)-iw*ih; + ov = iw*ih/ua; + } + return ov; +} diff --git a/random/include/Objectness/FilterTIG.h b/random/include/Objectness/FilterTIG.h new file mode 100755 index 0000000..101902c --- /dev/null +++ b/random/include/Objectness/FilterTIG.h @@ -0,0 +1,45 @@ +#pragma once + +class FilterTIG +{ +public: + void update(CMat &w); + + // For a W by H gradient magnitude map, find a W-7 by H-7 CV_32F matching score map + Mat matchTemplate(const Mat &mag1u); + + inline float dot(const INT64 tig1, const INT64 tig2, const INT64 tig4, const INT64 tig8); + +public: + void reconstruct(Mat &w); // For illustration purpose + +private: + static const int NUM_COMP = 2; // Number of components + static const int D = 64; // Dimension of TIG + INT64 _bTIGs[NUM_COMP]; // Binary TIG features + float _coeffs1[NUM_COMP]; // Coefficients of binary TIG features + + // For efficiently deals with different bits in CV_8U gradient map + float _coeffs2[NUM_COMP], _coeffs4[NUM_COMP], _coeffs8[NUM_COMP]; +}; + + +inline float FilterTIG::dot(const INT64 tig1, const INT64 tig2, const INT64 tig4, const INT64 tig8) +{ + INT64 bcT1 = __popcnt64(tig1); + INT64 bcT2 = __popcnt64(tig2); + INT64 bcT4 = __popcnt64(tig4); + INT64 bcT8 = __popcnt64(tig8); + + INT64 bc01 = (__popcnt64(_bTIGs[0] & tig1) << 1) - bcT1; + INT64 bc02 = ((__popcnt64(_bTIGs[0] & tig2) << 1) - bcT2) << 1; + INT64 bc04 = ((__popcnt64(_bTIGs[0] & tig4) << 1) - bcT4) << 2; + INT64 bc08 = ((__popcnt64(_bTIGs[0] & tig8) << 1) - bcT8) << 3; + + INT64 bc11 = (__popcnt64(_bTIGs[1] & tig1) << 1) - bcT1; + INT64 bc12 = ((__popcnt64(_bTIGs[1] & tig2) << 1) - bcT2) << 1; + INT64 bc14 = ((__popcnt64(_bTIGs[1] & tig4) << 1) - bcT4) << 2; + INT64 bc18 = ((__popcnt64(_bTIGs[1] & tig8) << 1) - bcT8) << 3; + + return _coeffs1[0] * (bc01 + bc02 + bc04 + bc08) + _coeffs1[1] * (bc11 + bc12 + bc14 + bc18); +} diff --git a/random/include/Objectness/ImgContrastBB.h b/random/include/Objectness/ImgContrastBB.h new file mode 100755 index 0000000..71c5fb5 --- /dev/null +++ b/random/include/Objectness/ImgContrastBB.h @@ -0,0 +1,62 @@ +#pragma once + +struct ImgContrastBB +{ + ImgContrastBB(CStr &imgPath); + ImgContrastBB(CMat &img3u); + + inline float contrastVal(Vec4i ¢er); + inline int regSum(Vec4i &box, Vec3i &sumColor); // Return region size and sum color + +private: + Mat iImg; + int _w, _h; + inline void assertBBox(Vec4i ¢er, CStr &name); +}; + +ImgContrastBB::ImgContrastBB(CStr &imgPath) +{ + Mat img3u = imread(imgPath); + integral(img3u, iImg, CV_32SC3); + _w = img3u.cols; + _h = img3u.rows; +} + +ImgContrastBB::ImgContrastBB(CMat &img3u) +{ + integral(img3u, iImg, CV_32SC3); + _w = img3u.cols; + _h = img3u.rows; +} + +int ImgContrastBB::regSum(Vec4i &box, Vec3i &sumColor) +{ + int x1 = box[0] - 1, y1 = box[1] - 1, x2 = box[2] - 1, y2 = box[3] - 1; + sumColor = iImg.at(y2, x2) + iImg.at(y1, x1) - iImg.at(y1, x2) - iImg.at(y2, x1); + return (x2 - x1)*(y2 - y1); +} + + +float ImgContrastBB::contrastVal(Vec4i ¢er) +{ + int wd = (center[2] - center[0])/2, hd = (center[3] - center[1])/2; + Vec4i surround(max(center[0] - wd, 1), max(center[1] - hd, 1), min(center[2] + wd, _w), min(center[3] + hd, _h)); + Vec3i cColr, sColr; + + assertBBox(center, "Center"); + assertBBox(center, "Surround"); + int cSz = regSum(center, cColr); + int sSz = regSum(surround, sColr); + + sColr -= cColr; + sSz -= cSz; + sColr /= sSz; + cColr /= cSz; + return sqrtf((float)(sqr(sColr[0] - cColr[0]) + sqr(sColr[1] - cColr[1]) + sqr(sColr[2] - cColr[2])))/100.0f; +} + +void ImgContrastBB::assertBBox(Vec4i ¢er, CStr &name) +{ + if (center[0] < 1 || center[1] < 1 || center[2] > _w || center[3] > _h) + printf("%s: (%d, %d, %d, %d), (%d, %d)\n", _S(name), center[0], center[1], center[2], center[3], _w, _h); +} \ No newline at end of file diff --git a/random/include/Objectness/Objectness.h b/random/include/Objectness/Objectness.h new file mode 100755 index 0000000..cc13fb8 --- /dev/null +++ b/random/include/Objectness/Objectness.h @@ -0,0 +1,74 @@ +#pragma once +#include "DataSetVOC.h" +#include "ValStructVec.h" +#include "FilterTIG.h" + +class Objectness +{ +public: + // base for window size quantization, feature window size (W, W), and non-maximal suppress size NSS + Objectness(double base = 2, int W = 8, int NSS = 2); + ~Objectness(void); + + // Load trained model. + int loadTrainedModel(string modelName = ""); // Return -1, 0, or 1 if partial, none, or all loaded + + // Get potential bounding boxes, each of which is represented by a Vec4i for (minX, minY, maxX, maxY). + // The trained model should be prepared before calling this function: loadTrainedModel() or trainStageI() + trainStageII(). + // Use numDet to control the final number of proposed bounding boxes, and number of per size (scale and aspect ratio) + void getObjBndBoxes(CMat &img3u, ValStructVec &valBoxes, int numDetPerSize = 120); + + // Read matrix from binary file + static bool matRead( const string& filename, Mat& M); + + static Mat aFilter(float delta, int sz); + +private: // Parameters + const double _base, _logBase; // base for window size quantization + const int _W; // As described in the paper: #Size, Size(_W, _H) of feature window. + const int _NSS; // Size for non-maximal suppress + const int _maxT, _minT, _numT; // The minimal and maximal dimensions of the template + + vecI _svmSzIdxs; // Indexes of active size. It's equal to _svmFilters.size() and _svmReW1f.rows + Mat _svmFilter; // Filters learned at stage I, each is a _H by _W CV_32F matrix + FilterTIG _tigF; // TIG filter + Mat _svmReW1f; // Re-weight parameters learned at stage II. + +private: // Help functions + + bool filtersLoaded() {int n = _svmSzIdxs.size(); return n > 0 && _svmReW1f.size() == Size(2, n) && _svmFilter.size() == Size(_W, _W);} + + int gtBndBoxSampling(const Vec4i &bbgt, vector &samples, vecI &bbR); + + Mat getFeature(CMat &img3u, const Vec4i &bb); // Return region feature + + inline double maxIntUnion(const Vec4i &bb, const vector &bbgts) {double maxV = 0; for(size_t i = 0; i < bbgts.size(); i++) maxV = max(maxV, DataSetVOC::interUnio(bb, bbgts[i])); return maxV; } + + // Convert VOC bounding box type to OpenCV Rect + inline Rect pnt2Rect(const Vec4i &bb){int x = bb[0] - 1, y = bb[1] - 1; return Rect(x, y, bb[2] - x, bb[3] - y);} + + // Template length at quantized scale t + inline int tLen(int t){return cvRound(pow(_base, t));} + + // Sub to quantization index + inline int sz2idx(int w, int h) {w -= _minT; h -= _minT; CV_Assert(w >= 0 && h >= 0 && w < _numT && h < _numT); return h * _numT + w + 1; } + inline string strVec4i(const Vec4i &v) const {return format("%d, %d, %d, %d", v[0], v[1], v[2], v[3]);} + + void predictBBoxSI(CMat &mag3u, ValStructVec &valBoxes, vecI &sz, int NUM_WIN_PSZ = 100, bool fast = true); + void predictBBoxSII(ValStructVec &valBoxes, const vecI &sz); + + // Calculate the image gradient: center option as in VLFeat + void gradientMag(CMat &imgBGR3u, Mat &mag1u); + + static void gradientRGB(CMat &bgr3u, Mat &mag1u); + static void gradientGray(CMat &bgr3u, Mat &mag1u); + static void gradientHSV(CMat &bgr3u, Mat &mag1u); + static void gradientXY(CMat &x1i, CMat &y1i, Mat &mag1u); + + static inline int bgrMaxDist(const Vec3b &u, const Vec3b &v) {int b = abs(u[0]-v[0]), g = abs(u[1]-v[1]), r = abs(u[2]-v[2]); b = max(b,g); return max(b,r);} + static inline int vecDist3b(const Vec3b &u, const Vec3b &v) {return abs(u[0]-v[0]) + abs(u[1]-v[1]) + abs(u[2]-v[2]);} + + //Non-maximal suppress + static void nonMaxSup(CMat &matchCost1f, ValStructVec &matchCost, int NSS = 1, int maxPoint = 50, bool fast = true); +}; + diff --git a/random/include/Objectness/ValStructVec.h b/random/include/Objectness/ValStructVec.h new file mode 100755 index 0000000..6a83a93 --- /dev/null +++ b/random/include/Objectness/ValStructVec.h @@ -0,0 +1,72 @@ +#pragma once + +/************************************************************************/ +/* A value struct vector that supports efficient sorting */ +/************************************************************************/ + +template +struct ValStructVec +{ + ValStructVec(){clear();} + inline int size() const {return sz;} + inline void clear() {sz = 0; structVals.clear(); valIdxes.clear();} + inline void reserve(int resSz){clear(); structVals.reserve(resSz); valIdxes.reserve(resSz); } + inline void pushBack(const VT& val, const ST& structVal) {valIdxes.push_back(make_pair(val, sz)); structVals.push_back(structVal); sz++;} + + inline const VT& operator ()(int i) const {return valIdxes[i].first;} // Should be called after sort + inline const ST& operator [](int i) const {return structVals[valIdxes[i].second];} // Should be called after sort + inline VT& operator ()(int i) {return valIdxes[i].first;} // Should be called after sort + inline ST& operator [](int i) {return structVals[valIdxes[i].second];} // Should be called after sort + + void sort(bool descendOrder = true); + const vector &getSortedStructVal(); + void append(const ValStructVec &newVals, int startV = 0); + + vector structVals; // struct values + +private: + int sz; // size of the value struct vector + vector > valIdxes; // Indexes after sort + bool smaller() {return true;}; + vector sortedStructVals; +}; + +template +void ValStructVec::append(const ValStructVec &newVals, int startV) +{ + int sz = newVals.size(); + for (int i = 0; i < sz; i++) + pushBack((float)((i+300)*startV)/*newVals(i)*/, newVals[i]); +} + +template +void ValStructVec::sort(bool descendOrder /* = true */) +{ + if (descendOrder) + std::sort(valIdxes.begin(), valIdxes.end(), std::greater >()); + else + std::sort(valIdxes.begin(), valIdxes.end(), std::less >()); +} + +template +const vector& ValStructVec::getSortedStructVal() +{ + sortedStructVals.resize(sz); + for (int i = 0; i < sz; i++) + sortedStructVals[i] = structVals[valIdxes[i].second]; + return sortedStructVals; +} + +/* +void valStructVecDemo() +{ + ValStructVec sVals; + sVals.pushBack(3, "String 3"); + sVals.pushBack(5, "String 5"); + sVals.pushBack(4, "String 4"); + sVals.pushBack(1, "String 1"); + sVals.sort(false); + for (int i = 0; i < sVals.size(); i++) + printf("%d, %s\n", sVals(i), _S(sVals[i])); +} +*/ diff --git a/random/include/Objectness/stdafx.h b/random/include/Objectness/stdafx.h new file mode 100755 index 0000000..a232b13 --- /dev/null +++ b/random/include/Objectness/stdafx.h @@ -0,0 +1,108 @@ +// stdafx.h : include file for standard system include files, +// or project specific include files that are used frequently, but +// are changed infrequently +// + +#pragma once + +#define _CRTDBG_MAP_ALLOC +#include +#ifdef _WIN32 +#include + +#include +#endif +#include + +#pragma once +#pragma warning(disable: 4996) +#pragma warning(disable: 4995) +#pragma warning(disable: 4805) +#pragma warning(disable: 4267) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +//#include + +#ifdef _WIN32 +#include +#include +#endif + +#ifndef __APPLE__ +#include +#endif +//#include + +#ifndef _WIN32 +#include +#include +#endif + +using namespace std; + +// TODO: reference additional headers your program requires here +#include "LibLinear/linear.h" +#include + +#define CV_VERSION_ID CVAUX_STR(CV_MAJOR_VERSION) CVAUX_STR(CV_MINOR_VERSION) CVAUX_STR(CV_SUBMINOR_VERSION) +#ifdef _DEBUG +#define cvLIB(name) "opencv_" name CV_VERSION_ID "d" +#else +#define cvLIB(name) "opencv_" name CV_VERSION_ID +#endif + +#pragma comment( lib, cvLIB("core")) +#pragma comment( lib, cvLIB("imgproc")) +#pragma comment( lib, cvLIB("highgui")) +using namespace cv; + +typedef vector vecI; +typedef const string CStr; +typedef const Mat CMat; +typedef vector vecS; +typedef vector vecM; +typedef vector vecF; +typedef vector vecD; + +#ifndef _WIN32 +typedef uint64_t UINT64; +typedef uint32_t UINT; +typedef int64_t INT64; +typedef uint8_t byte; +typedef bool BOOL; + +#define FALSE false; +#define __popcnt64 __builtin_popcountll +#endif + +enum{CV_FLIP_BOTH = -1, CV_FLIP_VERTICAL = 0, CV_FLIP_HORIZONTAL = 1}; +#define _S(str) ((str).c_str()) +#define CHK_IND(p) ((p).x >= 0 && (p).x < _w && (p).y >= 0 && (p).y < _h) +#define CV_Assert_(expr, args) \ +{\ + if(!(expr)) {\ + string msg = cv::format args; \ + printf("%s in %s:%d\n", msg.c_str(), __FILE__, __LINE__); \ + cv::error(cv::Exception(CV_StsAssert, msg, __FUNCTION__, __FILE__, __LINE__) ); }\ +} + +// Return -1 if not in the list +template +static inline int findFromList(const T &word, const vector &strList) {size_t idx = find(strList.begin(), strList.end(), word) - strList.begin(); return idx < strList.size() ? idx : -1;} +template inline T sqr(T x) { return x * x; } // out of range risk for T = byte, ... +template inline T vecSqrDist(const Vec &v1, const Vec &v2) {T s = 0; for (int i=0; i inline T vecDist(const Vec &v1, const Vec &v2) { return sqrt(vecSqrDist(v1, v2)); } // out of range risk for T = byte, ... + +inline Rect Vec4i2Rect(Vec4i &v){return Rect(Point(v[0] - 1, v[1] - 1), Point(v[2], v[3])); } + +#include "CmFile.h" +#include "CmTimer.h" + diff --git a/random/lib/libobjectness.dylib b/random/lib/libobjectness.dylib new file mode 100755 index 0000000..cc66810 Binary files /dev/null and b/random/lib/libobjectness.dylib differ diff --git a/random/model/ObjNessB2W8MAXBGR.idx b/random/model/ObjNessB2W8MAXBGR.idx new file mode 100644 index 0000000..1b8e988 Binary files /dev/null and b/random/model/ObjNessB2W8MAXBGR.idx differ diff --git a/random/model/ObjNessB2W8MAXBGR.wS1 b/random/model/ObjNessB2W8MAXBGR.wS1 new file mode 100644 index 0000000..88b0ef0 Binary files /dev/null and b/random/model/ObjNessB2W8MAXBGR.wS1 differ diff --git a/random/model/ObjNessB2W8MAXBGR.wS2 b/random/model/ObjNessB2W8MAXBGR.wS2 new file mode 100644 index 0000000..3fa5e32 Binary files /dev/null and b/random/model/ObjNessB2W8MAXBGR.wS2 differ diff --git a/random/model/ObjNessB2W8MAXBGR.xN b/random/model/ObjNessB2W8MAXBGR.xN new file mode 100644 index 0000000..466ecf1 Binary files /dev/null and b/random/model/ObjNessB2W8MAXBGR.xN differ diff --git a/random/model/ObjNessB2W8MAXBGR.xP b/random/model/ObjNessB2W8MAXBGR.xP new file mode 100644 index 0000000..01bdefd Binary files /dev/null and b/random/model/ObjNessB2W8MAXBGR.xP differ diff --git a/random/src/.main.cpp.swo b/random/src/.main.cpp.swo new file mode 100644 index 0000000..ef40996 Binary files /dev/null and b/random/src/.main.cpp.swo differ diff --git a/random/src/CVBoostConverter.hpp b/random/src/CVBoostConverter.hpp new file mode 100644 index 0000000..5bbf456 --- /dev/null +++ b/random/src/CVBoostConverter.hpp @@ -0,0 +1,423 @@ +/* + * CVBoostConverter.hpp + * + * Created on: Mar 20, 2014 + * Author: Gregory Kramida + * Copyright: (c) 2014 Gregory Kramida + * License: MIT + */ + +#ifndef CVBOOSTCONVERTER_HPP_ +#define CVBOOSTCONVERTER_HPP_ + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include +#include +#include +#include +#include +#include + +using namespace cv; +namespace bcvt{ +static PyObject* opencv_error = 0; + +class PyAllowThreads { +public: + PyAllowThreads() : + _state(PyEval_SaveThread()) { + } + ~PyAllowThreads() { + PyEval_RestoreThread(_state); + } +private: + PyThreadState* _state; +}; + +class PyEnsureGIL { +public: + PyEnsureGIL() : + _state(PyGILState_Ensure()) { + } + ~PyEnsureGIL() { + PyGILState_Release(_state); + } +private: + PyGILState_STATE _state; +}; +#define ERRWRAP2(expr) \ +try \ +{ \ + PyAllowThreads allowThreads; \ + expr; \ +} \ +catch (const cv::Exception &e) \ +{ \ + PyErr_SetString(opencv_error, e.what()); \ + return 0; \ +} + +class NumpyAllocator: public MatAllocator { +public: + NumpyAllocator() { + stdAllocator = Mat::getStdAllocator(); + } + ~NumpyAllocator() { + } + + UMatData* allocate(PyObject* o, int dims, const int* sizes, int type, + size_t* step) const { + UMatData* u = new UMatData(this); + u->data = u->origdata = (uchar*) PyArray_DATA((PyArrayObject*) o); + npy_intp* _strides = PyArray_STRIDES((PyArrayObject*) o); + for (int i = 0; i < dims - 1; i++) + step[i] = (size_t) _strides[i]; + step[dims - 1] = CV_ELEM_SIZE(type); + u->size = sizes[0] * step[0]; + u->userdata = o; + return u; + } + + UMatData* allocate(int dims0, const int* sizes, int type, void* data, + size_t* step, int flags, UMatUsageFlags usageFlags) const { + if (data != 0) { + CV_Error(Error::StsAssert, "The data should normally be NULL!"); + // probably this is safe to do in such extreme case + return stdAllocator->allocate(dims0, sizes, type, data, step, flags, + usageFlags); + } + PyEnsureGIL gil; + + int depth = CV_MAT_DEPTH(type); + int cn = CV_MAT_CN(type); + const int f = (int) (sizeof(size_t) / 8); + int typenum = + depth == CV_8U ? NPY_UBYTE : + depth == CV_8S ? NPY_BYTE : + depth == CV_16U ? NPY_USHORT : + depth == CV_16S ? NPY_SHORT : + depth == CV_32S ? NPY_INT : + depth == CV_32F ? NPY_FLOAT : + depth == CV_64F ? + NPY_DOUBLE : f * NPY_ULONGLONG + (f ^ 1) * NPY_UINT; + int i, dims = dims0; + cv::AutoBuffer _sizes(dims + 1); + for (i = 0; i < dims; i++) + _sizes[i] = sizes[i]; + if (cn > 1) + _sizes[dims++] = cn; + PyObject* o = PyArray_SimpleNew(dims, _sizes, typenum); + if (!o) + CV_Error_(Error::StsError, + ("The numpy array of typenum=%d, ndims=%d can not be created", typenum, dims)); + return allocate(o, dims0, sizes, type, step); + } + + bool allocate(UMatData* u, int accessFlags, + UMatUsageFlags usageFlags) const { + return stdAllocator->allocate(u, accessFlags, usageFlags); + } + + void deallocate(UMatData* u) const { + if (u) { + PyEnsureGIL gil; + PyObject* o = (PyObject*) u->userdata; + Py_XDECREF(o); + delete u; + } + } + + const MatAllocator* stdAllocator; +}; + +NumpyAllocator g_numpyAllocator; + +PyObject* fromMatToNDArray(const Mat& m) { + if (!m.data) + Py_RETURN_NONE; + Mat temp, + *p = (Mat*) &m; + if (!p->u || p->allocator != &g_numpyAllocator) { + temp.allocator = &g_numpyAllocator; + ERRWRAP2(m.copyTo(temp)); + p = &temp; + } + PyObject* o = (PyObject*) p->u->userdata; + Py_INCREF(o); + return o; +} + +static int failmsg(const char *fmt, ...) { + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +Mat fromNDArrayToMat(PyObject* o) { + cv::Mat m; + bool allowND = true; + if (!PyArray_Check(o)) { + failmsg("argument is not a numpy array"); + if (!m.data) + m.allocator = &g_numpyAllocator; + } else { + PyArrayObject* oarr = (PyArrayObject*) o; + + bool needcopy = false, needcast = false; + int typenum = PyArray_TYPE(oarr), new_typenum = typenum; + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : + typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT ? CV_32S : + typenum == NPY_INT32 ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if (type < 0) { + if (typenum == NPY_INT64 || typenum == NPY_UINT64 + || type == NPY_LONG) { + needcopy = needcast = true; + new_typenum = NPY_INT; + type = CV_32S; + } else { + failmsg("Argument data type is not supported"); + m.allocator = &g_numpyAllocator; + return m; + } + } + +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + + int ndims = PyArray_NDIM(oarr); + if (ndims >= CV_MAX_DIM) { + failmsg("Dimensionality of argument is too high"); + if (!m.data) + m.allocator = &g_numpyAllocator; + return m; + } + + int size[CV_MAX_DIM + 1]; + size_t step[CV_MAX_DIM + 1]; + size_t elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(oarr); + const npy_intp* _strides = PyArray_STRIDES(oarr); + bool ismultichannel = ndims == 3 && _sizes[2] <= CV_CN_MAX; + + for (int i = ndims - 1; i >= 0 && !needcopy; i--) { + // these checks handle cases of + // a) multi-dimensional (ndims > 2) arrays, as well as simpler 1- and 2-dimensional cases + // b) transposed arrays, where _strides[] elements go in non-descending order + // c) flipped arrays, where some of _strides[] elements are negative + if ((i == ndims - 1 && (size_t) _strides[i] != elemsize) + || (i < ndims - 1 && _strides[i] < _strides[i + 1])) + needcopy = true; + } + + if (ismultichannel && _strides[1] != (npy_intp) elemsize * _sizes[2]) + needcopy = true; + + if (needcopy) { + + if (needcast) { + o = PyArray_Cast(oarr, new_typenum); + oarr = (PyArrayObject*) o; + } else { + oarr = PyArray_GETCONTIGUOUS(oarr); + o = (PyObject*) oarr; + } + + _strides = PyArray_STRIDES(oarr); + } + + for (int i = 0; i < ndims; i++) { + size[i] = (int) _sizes[i]; + step[i] = (size_t) _strides[i]; + } + + // handle degenerate case + if (ndims == 0) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if (ismultichannel) { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + + if (ndims > 2 && !allowND) { + failmsg("%s has more than 2 dimensions"); + } else { + + m = Mat(ndims, size, type, PyArray_DATA(oarr), step); + m.u = g_numpyAllocator.allocate(o, ndims, size, type, step); + m.addref(); + + if (!needcopy) { + Py_INCREF(o); + } + } + m.allocator = &g_numpyAllocator; + } + return m; +} + +struct matToNDArrayBoostConverter { + static PyObject* convert(Mat const& m) { + if (!m.data) + Py_RETURN_NONE; + Mat *p = (Mat*) &m; + Mat temp; + if (!p->u || p->allocator != &g_numpyAllocator) { + temp.allocator = &g_numpyAllocator; + ERRWRAP2(m.copyTo(temp)); + p = &temp; + } + PyObject* o = (PyObject*) p->u->userdata; + return boost::python::incref(o); + } +}; + +struct matFromNDArrayBoostConverter { + + matFromNDArrayBoostConverter() { + boost::python::converter::registry::push_back(convertible, construct, + boost::python::type_id()); + } + + /// @brief Check if PyObject is an array and can be converted to OpenCV matrix. + static void* convertible(PyObject* object) { + if (!PyArray_Check(object)) { + return NULL; + } +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + PyArrayObject* oarr = (PyArrayObject*) object; + + int typenum = PyArray_TYPE(oarr); + if (typenum != NPY_INT64 && typenum != NPY_UINT64 && typenum != NPY_LONG + && typenum != NPY_UBYTE && typenum != NPY_BYTE + && typenum != NPY_USHORT && typenum != NPY_SHORT + && typenum != NPY_INT && typenum != NPY_INT32 + && typenum != NPY_FLOAT && typenum != NPY_DOUBLE) { + return NULL; + } + int ndims = PyArray_NDIM(oarr); //data type not supported + + if (ndims >= CV_MAX_DIM) { + return NULL; //too many dimensions + } + return object; + } + + /// @brief Construct a Mat from an NDArray object. + static void construct(PyObject* object, + boost::python::converter::rvalue_from_python_stage1_data* data) { + namespace python = boost::python; + // Object is a borrowed reference, so create a handle indicting it is + // borrowed for proper reference counting. + python::handle<> handle(python::borrowed(object)); + + // Obtain a handle to the memory block that the converter has allocated + // for the C++ type. + typedef python::converter::rvalue_from_python_storage storage_type; + void* storage = reinterpret_cast(data)->storage.bytes; + + // Allocate the C++ type into the converter's memory block, and assign + // its handle to the converter's convertible variable. The C++ + // container is populated by passing the begin and end iterators of + // the python object to the container's constructor. + PyArrayObject* oarr = (PyArrayObject*) object; + + bool needcopy = false, needcast = false; + int typenum = PyArray_TYPE(oarr), new_typenum = typenum; + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : + typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT ? CV_32S : + typenum == NPY_INT32 ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if (type < 0) { + needcopy = needcast = true; + new_typenum = NPY_INT; + type = CV_32S; + } + +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + int ndims = PyArray_NDIM(oarr); + + int size[CV_MAX_DIM + 1]; + size_t step[CV_MAX_DIM + 1]; + size_t elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(oarr); + const npy_intp* _strides = PyArray_STRIDES(oarr); + bool ismultichannel = ndims == 3 && _sizes[2] <= CV_CN_MAX; + + for (int i = ndims - 1; i >= 0 && !needcopy; i--) { + // these checks handle cases of + // a) multi-dimensional (ndims > 2) arrays, as well as simpler 1- and 2-dimensional cases + // b) transposed arrays, where _strides[] elements go in non-descending order + // c) flipped arrays, where some of _strides[] elements are negative + if ((i == ndims - 1 && (size_t) _strides[i] != elemsize) + || (i < ndims - 1 && _strides[i] < _strides[i + 1])) + needcopy = true; + } + + if (ismultichannel && _strides[1] != (npy_intp) elemsize * _sizes[2]) + needcopy = true; + + if (needcopy) { + + if (needcast) { + object = PyArray_Cast(oarr, new_typenum); + oarr = (PyArrayObject*) object; + } else { + oarr = PyArray_GETCONTIGUOUS(oarr); + object = (PyObject*) oarr; + } + + _strides = PyArray_STRIDES(oarr); + } + + for (int i = 0; i < ndims; i++) { + size[i] = (int) _sizes[i]; + step[i] = (size_t) _strides[i]; + } + + // handle degenerate case + if (ndims == 0) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if (ismultichannel) { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + if (!needcopy) { + Py_INCREF(object); + } + cv::Mat* m = new (storage) cv::Mat(ndims, size, type, PyArray_DATA(oarr), step); + m->u = g_numpyAllocator.allocate(object, ndims, size, type, step); + m->allocator = &g_numpyAllocator; + m->addref(); + data->convertible = storage; + } +}; +} // end namespace bcvt +#endif /* CVBOOSTCONVERTER_HPP_ */ diff --git a/random/src/adjacency.h b/random/src/adjacency.h new file mode 100644 index 0000000..dd4e3ee --- /dev/null +++ b/random/src/adjacency.h @@ -0,0 +1,26 @@ +#pragma once + +#include + +class AdjacencyMatrix : public std::vector { +public: + AdjacencyMatrix(int n) { + this->n = n; + assign((n + 1)*n/2, false); + } + + reference get(int i, int j) { + return at(index(i, j)); + } + +protected: + int n; + + int index(int i, int j) { + int tmp = i > j ? j : i; + j = i > j ? i : j; + i = tmp; + return (n * i) + j - ((i * (i+1)) / 2); + } +}; + diff --git a/random/src/connection.h b/random/src/connection.h new file mode 100644 index 0000000..ccd48a4 --- /dev/null +++ b/random/src/connection.h @@ -0,0 +1,25 @@ +#pragma once + +#include "segment.h" + +class Connection +{ +public: + Connection(const Segment & a_, const Segment & b_, uint8_t flags) : + a(a_.id), b(b_.id), similarity(exp(Segment::computeSimilarity(a_, b_, flags))) + { + } + + friend std::ostream& operator<<(std::ostream &out, Connection & c) { + return out << ""; + } + + bool operator<(const Connection & rhs) const { + return rhs.similarity < similarity; + } + + int a; + int b; + float similarity; +}; + diff --git a/random/src/conversion.cpp b/random/src/conversion.cpp new file mode 100644 index 0000000..475c670 --- /dev/null +++ b/random/src/conversion.cpp @@ -0,0 +1,231 @@ +# include "conversion.h" +/* + * The following conversion functions are taken/adapted from OpenCV's cv2.cpp file + * inside modules/python/src2 folder. + */ + +static void init() +{ + import_array(); +} + +static int failmsg(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class PyAllowThreads +{ +public: + PyAllowThreads() : _state(PyEval_SaveThread()) {} + ~PyAllowThreads() + { + PyEval_RestoreThread(_state); + } +private: + PyThreadState* _state; +}; + +class PyEnsureGIL +{ +public: + PyEnsureGIL() : _state(PyGILState_Ensure()) {} + ~PyEnsureGIL() + { + PyGILState_Release(_state); + } +private: + PyGILState_STATE _state; +}; + +using namespace cv; + +static PyObject* failmsgp(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class NumpyAllocator : public MatAllocator +{ +public: + NumpyAllocator() {} + ~NumpyAllocator() {} + + void allocate(int dims, const int* sizes, int type, int*& refcount, + uchar*& datastart, uchar*& data, size_t* step) + { + PyEnsureGIL gil; + + int depth = CV_MAT_DEPTH(type); + int cn = CV_MAT_CN(type); + const int f = (int)(sizeof(size_t)/8); + int typenum = depth == CV_8U ? NPY_UBYTE : depth == CV_8S ? NPY_BYTE : + depth == CV_16U ? NPY_USHORT : depth == CV_16S ? NPY_SHORT : + depth == CV_32S ? NPY_INT : depth == CV_32F ? NPY_FLOAT : + depth == CV_64F ? NPY_DOUBLE : f*NPY_ULONGLONG + (f^1)*NPY_UINT; + int i; + npy_intp _sizes[CV_MAX_DIM+1]; + for( i = 0; i < dims; i++ ) + { + _sizes[i] = sizes[i]; + } + + if( cn > 1 ) + { + _sizes[dims++] = cn; + } + + PyObject* o = PyArray_SimpleNew(dims, _sizes, typenum); + + if(!o) + { + CV_Error_(CV_StsError, ("The numpy array of typenum=%d, ndims=%d can not be created", typenum, dims)); + } + refcount = refcountFromPyObject(o); + + npy_intp* _strides = PyArray_STRIDES(o); + for( i = 0; i < dims - (cn > 1); i++ ) + step[i] = (size_t)_strides[i]; + datastart = data = (uchar*)PyArray_DATA(o); + } + + void deallocate(int* refcount, uchar*, uchar*) + { + PyEnsureGIL gil; + if( !refcount ) + return; + PyObject* o = pyObjectFromRefcount(refcount); + Py_INCREF(o); + Py_DECREF(o); + } +}; + +NumpyAllocator g_numpyAllocator; + +NDArrayConverter::NDArrayConverter() { init(); } + +void NDArrayConverter::init() +{ + import_array(); +} + +cv::Mat NDArrayConverter::toMat(const PyObject *o) +{ + cv::Mat m; + + if(!o || o == Py_None) + { + if( !m.data ) + m.allocator = &g_numpyAllocator; + } + + if( !PyArray_Check(o) ) + { + failmsg("toMat: Object is not a numpy array"); + } + + int typenum = PyArray_TYPE(o); + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT || typenum == NPY_LONG ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if( type < 0 ) + { + failmsg("toMat: Data type = %d is not supported", typenum); + } + + int ndims = PyArray_NDIM(o); + + if(ndims >= CV_MAX_DIM) + { + failmsg("toMat: Dimensionality (=%d) is too high", ndims); + } + + int size[CV_MAX_DIM+1]; + size_t step[CV_MAX_DIM+1], elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(o); + const npy_intp* _strides = PyArray_STRIDES(o); + bool transposed = false; + + for(int i = 0; i < ndims; i++) + { + size[i] = (int)_sizes[i]; + step[i] = (size_t)_strides[i]; + } + + if( ndims == 0 || step[ndims-1] > elemsize ) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if( ndims >= 2 && step[0] < step[1] ) + { + std::swap(size[0], size[1]); + std::swap(step[0], step[1]); + transposed = true; + } + + if( ndims == 3 && size[2] <= CV_CN_MAX && step[1] == elemsize*size[2] ) + { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + + if( ndims > 2) + { + failmsg("toMat: Object has more than 2 dimensions"); + } + + m = Mat(ndims, size, type, PyArray_DATA(o), step); + + if( m.data ) + { + m.refcount = refcountFromPyObject(o); + m.addref(); // protect the original numpy array from deallocation + // (since Mat destructor will decrement the reference counter) + }; + m.allocator = &g_numpyAllocator; + + if( transposed ) + { + Mat tmp; + tmp.allocator = &g_numpyAllocator; + transpose(m, tmp); + m = tmp; + } + return m; +} + +PyObject* NDArrayConverter::toNDArray(const cv::Mat& m) +{ + if( !m.data ) + Py_RETURN_NONE; + Mat temp, *p = (Mat*)&m; + if(!p->refcount || p->allocator != &g_numpyAllocator) + { + temp.allocator = &g_numpyAllocator; + m.copyTo(temp); + p = &temp; + } + p->addref(); + return pyObjectFromRefcount(p->refcount); +} diff --git a/random/src/conversion.h b/random/src/conversion.h new file mode 100644 index 0000000..3975e17 --- /dev/null +++ b/random/src/conversion.h @@ -0,0 +1,60 @@ +# ifndef __COVERSION_OPENCV_H__ +# define __COVERSION_OPENCV_H__ + +#include +#include +#include +#include +#include "numpy/ndarrayobject.h" + +static PyObject* opencv_error = 0; + +static int failmsg(const char *fmt, ...); + +class PyAllowThreads; + +class PyEnsureGIL; + +#define ERRWRAP2(expr) \ +try \ +{ \ + PyAllowThreads allowThreads; \ + expr; \ +} \ +catch (const cv::Exception &e) \ +{ \ + PyErr_SetString(opencv_error, e.what()); \ + return 0; \ +} + +static PyObject* failmsgp(const char *fmt, ...); + +static size_t REFCOUNT_OFFSET = (size_t)&(((PyObject*)0)->ob_refcnt) + + (0x12345678 != *(const size_t*)"\x78\x56\x34\x12\0\0\0\0\0")*sizeof(int); + +static inline PyObject* pyObjectFromRefcount(const int* refcount) +{ + return (PyObject*)((size_t)refcount - REFCOUNT_OFFSET); +} + +static inline int* refcountFromPyObject(const PyObject* obj) +{ + return (int*)((size_t)obj + REFCOUNT_OFFSET); +} + + +class NumpyAllocator; + +enum { ARG_NONE = 0, ARG_MAT = 1, ARG_SCALAR = 2 }; + +class NDArrayConverter +{ +private: + void init(); +public: + NDArrayConverter(); + cv::Mat toMat(const PyObject* o); + PyObject* toNDArray(const cv::Mat& mat); +}; + +# endif diff --git a/random/src/location_prior.h b/random/src/location_prior.h new file mode 100644 index 0000000..1587936 --- /dev/null +++ b/random/src/location_prior.h @@ -0,0 +1,44 @@ +#pragma once + +#include "selection_prior.h" + +const float deviation = 4.f; + +class LocationPrior : public SelectionPrior { +public: + virtual SelectionPriorMap computeSelectionPrior(const cv::Mat & image, const std::vector & segments) { + SelectionPriorMap prior; + + uint32_t radius = (image.cols*image.cols + image.rows*image.rows) / 4; + cv::Point center = cv::Point(image.cols / 2, image.rows / 2); + float sum = 0.f; + + for (const auto s: segments) { + if (s.empty()) + continue; + + cv::Point diff = (s.min_p + s.max_p) * 0.5 - center; + float likelihood = exp(-diff.dot(diff) / float(radius) * deviation); + sum += likelihood; + prior.insert({s.id, likelihood}); + } + + for (auto & p: prior) { + p.second /= sum; + } + + return prior; + +/* std::uniform_real_distribution dis(0.f, sum); + float rnd = dis(gen); + sum = 0.f; + for (auto p: prior) { + sum += p.second; + if (sum >= rnd) { + return p.first; + } + } + + return 0;*/ + } +}; diff --git a/random/src/main.cpp b/random/src/main.cpp new file mode 100644 index 0000000..e04f823 --- /dev/null +++ b/random/src/main.cpp @@ -0,0 +1,212 @@ +//#include "CVBoostConverter.hpp" +#include +#include "conversion.h" +#include + +#include "connection.h" +#include "adjacency.h" +#include "uniform_prior.h" +#include "location_prior.h" +#include "objectness_prior.h" +#include "objectness_location_prior.h" +#include "random_stopping_criterion.h" + +using namespace boost::python; + +//#define DEBUG 1 + +std::random_device rd_; +std::mt19937 gen_(rd_()); + +cv::Mat get_bboxes_(cv::Mat image, cv::Mat seg, cv::Mat edge, uint8_t flags, uint32_t n, std::string selection_prior) { + double max_id_; + cv::minMaxIdx(seg, nullptr, &max_id_); + int max_id = max_id_; + + std::vector segments; + segments.reserve(max_id); + int nchannels = image.channels(); + cv::Size size = image.size(); + for (int i = 0; i <= max_id; i++) { + segments.push_back(Segment(i, nchannels, size)); + } + + { + AdjacencyMatrix adjacency(max_id + 1); + for (int i = 0; i < image.rows; i++) { + for (int j = 0; j < image.cols; j++) { + cv::Point p(j, i); + uint16_t id = seg.at(p); + if (nchannels == 1) + segments[id].addPoint(image.at>(p), edge.at(p), p); + else if (nchannels == 3) + segments[id].addPoint(image.at(p), edge.at(p), p); + + if (i < image.rows - 1) { + uint16_t n = seg.at(i+1, j); + if (n != id && adjacency.get(id, n) == false) { + adjacency.get(id, n) = true; + segments[id].addNeighbour(n); + segments[n].addNeighbour(id); + } + } + + if (j < image.cols - 1) { + uint16_t n = seg.at(i, j+1); + if (n != id && adjacency.get(id, n) == false) { + adjacency.get(id, n) = true; + segments[id].addNeighbour(n); + segments[n].addNeighbour(id); + } + } + } + } + } + + cv::Mat bboxes; + float similarity_sum = 0.f; + for (auto & s: segments) { + s.normalizeHistogram(); + + cv::Mat bbox = cv::Mat(1, 4, CV_32SC1); + bbox.at(0) = s.min_p.x; + bbox.at(1) = s.min_p.y; + bbox.at(2) = s.max_p.x; + bbox.at(3) = s.max_p.y; + if (bboxes.empty()) + bboxes = bbox; + else + cv::vconcat(bboxes, bbox, bboxes); + } + + SelectionPriorMap prior; + + if (selection_prior == "uniform") { + UniformPrior up; + prior = up.computeSelectionPrior(image, segments); + } + else if (selection_prior == "location") { + LocationPrior lp; + prior = lp.computeSelectionPrior(image, segments); + } + else if (selection_prior == "objectness") { + ObjectnessPrior op; + prior = op.computeSelectionPrior(image, segments); + } + else if (selection_prior == "objectness-location") { + ObjectnessLocationPrior olp; + prior = olp.computeSelectionPrior(image, segments); + } + else { + throw std::runtime_error("Unknown selection prior method: " + selection_prior); + } + + Segment s; + for (uint32_t i = 0; i < n; i++) { + s = segments[prior.poll()]; + RandomStoppingCriterion stop; + cv::Rect r(s.min_p, s.max_p); + + while (stop.stop(image, r) == false) { + float sum = 0.f; + std::vector connections; + for (auto n: s.neighbours) { + connections.push_back(Connection(s, segments[n], flags)); + sum += (connections.end() - 1)->similarity; + } + + std::uniform_real_distribution dis(0.f, sum); + float rnd = dis(gen_); + + sum = 0.f; + for (auto & c: connections) { + sum += c.similarity; + if (sum >= rnd) { + s = Segment::merge(s, segments[c.b]); + break; + } + } + } + + cv::Mat bbox = cv::Mat(1, 4, CV_32SC1); + bbox.at(0) = s.min_p.x; + bbox.at(1) = s.min_p.y; + bbox.at(2) = s.max_p.x; + bbox.at(3) = s.max_p.y; + if (bboxes.empty()) + bboxes = bbox; + else + cv::vconcat(bboxes, bbox, bboxes); + } + +#ifdef DEBUG + cv::Mat draw; + image.copyTo(draw); + cv::cvtColor(draw, draw, cv::COLOR_BGR2GRAY); + cv::Rect r(cv::Point(bboxes.at(0), bboxes.at(1)), cv::Point(bboxes.at(2), bboxes.at(3))); + cv::rectangle(draw, r, cv::Scalar(255)); + cv::addWeighted(draw, 0.5, s.mask * 255, 0.5, 0.0, draw); + cv::imshow("Image", draw); + + std::cout << bboxes << std::endl; +#endif + + return bboxes; +} + +PyObject * get_bboxes(PyObject * image_, PyObject * seg_, PyObject * edge_, uint8_t flags, uint32_t n, std::string selection_prior) { + NDArrayConverter cvt; + cv::Mat image, seg, edge; + image = cvt.toMat(image_); + seg = cvt.toMat(seg_); + edge = cvt.toMat(edge_); + return cvt.toNDArray(get_bboxes_(image, seg, edge, flags, n, selection_prior)); +} + +static void init_ar() { + Py_Initialize(); + import_array(); +} + +BOOST_PYTHON_MODULE(random_selection) { + init_ar(); + + def("get_bboxes", get_bboxes); +} + +int main(int argc, char * argv[]) { + if (argc != 6) { + std::cout << "Usage: " << argv[0] << " " << std::endl; + return 0; + } + + cv::Mat image = cv::imread(argv[1]); + cv::Mat seg = cv::imread(argv[2], cv::IMREAD_UNCHANGED); + cv::Mat edge = cv::imread(argv[3], cv::IMREAD_UNCHANGED); + uint8_t iterations = atoi(argv[4]); + std::string selection_prior = std::string(argv[5]); + +// cv::namedWindow("Image", cv::WINDOW_NORMAL); +// cv::imshow("Image", (seg == 338) * 255); +// cv::waitKey(); + +// cv::namedWindow("Image", cv::WINDOW_NORMAL); +// while (cv::waitKey() != 'q') { + std::clock_t begin = std::clock(); + cv::Mat bboxes = get_bboxes_(image, seg, edge, COLOR_SIMILARITY | TEXTURE_SIMILARITY/* | SIZE_SIMILARITY | BBOX_SIMILARITY*/, iterations, "objectness"); + std::clock_t end = std::clock(); + double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC; + std::cout << "Times passed in seconds: " << elapsed_secs << std::endl; +// } + +// for (int i = 0; i < bboxes.rows; i++) { +// cv::Point p1(bboxes.at(i,0), bboxes.at(i,1)); +// cv::Point p2(bboxes.at(i,2), bboxes.at(i,3)); +// cv::Mat tmp; +// image.copyTo(tmp); +// cv::rectangle(tmp, p1, p2, cv::Scalar(255, 255, 0)); +// cv::imshow("Image", tmp); +// cv::waitKey(); +// } + return 0; +} diff --git a/random/src/objectness_location_prior.h b/random/src/objectness_location_prior.h new file mode 100644 index 0000000..8094cd4 --- /dev/null +++ b/random/src/objectness_location_prior.h @@ -0,0 +1,77 @@ +#pragma once + +#include "selection_prior.h" +#include "Objectness/stdafx.h" +#include "Objectness/Objectness.h" + +class ObjectnessLocationPrior : public SelectionPrior { +public: + ObjectnessLocationPrior() { + obj.loadTrainedModel("/Users/hans/MasterThesis/code/cpp/random/model/ObjNessB2W8MAXBGR"); + } + + virtual SelectionPriorMap computeSelectionPrior(const cv::Mat & image, const std::vector & segments) { + SelectionPriorMap prior; + + cv::Mat likelihood = cv::Mat::zeros(image.size(), CV_32FC1); + + ValStructVec boxes; + obj.getObjBndBoxes(image, boxes, 130); + +// cv::Mat test = cv::imread("../../../images/horse.jpg"); +// cv::Mat test; image.copyTo(test); +// for (int i = 0; i < boxes.size(); i++) { +// cv::Mat draw; +// test.copyTo(draw); +// cv::rectangle(draw, cv::Point(boxes[i][0], boxes[i][1]), cv::Point(boxes[i][2], boxes[i][3]), cv::Scalar(255, 255, 255)); +// std::cout << boxes(i) << std::endl; +// cv::namedWindow("Objectness", 0); +// cv::imshow("Objectness", draw); +// cv::waitKey(); +// } + + float max = boxes(0); + float min = boxes(boxes.size() - 1); + for (int i = 0; i < boxes.size(); i++) { + float score = (boxes(i) - min) / (max - min); // normalized score +// std::cout << "score: " << score << std::endl; + for (int x = boxes[i][0] - 1; x < boxes[i][2]; x++) { + for (int y = boxes[i][1] - 1; y < boxes[i][3]; y++) { + likelihood.at(y, x) += score; + } + } + } + + uint32_t radius = (image.cols*image.cols + image.rows*image.rows) / 4; + cv::Point center = cv::Point(image.cols / 2, image.rows / 2); + + float sum = 0.f; + for (const auto & s: segments) { + cv::Mat mask = cv::Mat::zeros(image.size(), CV_8UC1); + cv::rectangle(mask, s.min_p, s.max_p, cv::Scalar(1), CV_FILLED); + float score = cv::mean(likelihood, mask)[0]; + + cv::Point diff = (s.min_p + s.max_p) * 0.5 - center; + score *= exp(-diff.dot(diff) / float(radius) * deviation); + + sum += score; + prior.insert({s.id, score}); + } + + for (auto & p: prior) { + p.second /= sum; + } + +// cv::normalize(likelihood, likelihood, 0.f, 1.f, cv::NORM_MINMAX); +// cv::namedWindow("Likelihood", cv::WINDOW_AUTOSIZE); +// cv::namedWindow("Image", cv::WINDOW_AUTOSIZE); +// cv::imshow("Likelihood", likelihood); +// cv::imshow("Image", image); +// cv::waitKey(); + + return prior; + } + +protected: + Objectness obj; +}; diff --git a/random/src/objectness_prior.h b/random/src/objectness_prior.h new file mode 100644 index 0000000..f28e07d --- /dev/null +++ b/random/src/objectness_prior.h @@ -0,0 +1,85 @@ +#pragma once + +#include "selection_prior.h" +#include "Objectness/stdafx.h" +#include "Objectness/Objectness.h" + +class ObjectnessPrior : public SelectionPrior { +public: + ObjectnessPrior() { + obj.loadTrainedModel("/Users/hans/MasterThesis/code/cpp/random/model/ObjNessB2W8MAXBGR"); + } + + virtual SelectionPriorMap computeSelectionPrior(const cv::Mat & image, const std::vector & segments) { + SelectionPriorMap prior; + + cv::Mat likelihood = cv::Mat::zeros(image.size(), CV_32FC1); + + ValStructVec boxes; + obj.getObjBndBoxes(image, boxes, 130); + +// cv::Mat test = cv::imread("../../../images/horse.jpg"); +// cv::Mat test; image.copyTo(test); +// for (int i = 0; i < boxes.size(); i++) { +// cv::Mat draw; +// test.copyTo(draw); +// cv::rectangle(draw, cv::Point(boxes[i][0], boxes[i][1]), cv::Point(boxes[i][2], boxes[i][3]), cv::Scalar(255, 255, 255)); +// std::cout << boxes(i) << std::endl; +// cv::namedWindow("Objectness", 0); +// cv::imshow("Objectness", draw); +// cv::waitKey(); +// } + + float sum = 0.f; + for (const auto & s: segments) { + float score = 0.f; + for (int i = 0; i < boxes.size(); i++) { + const cv::Vec4i & box = boxes[i]; + if (s.min_p.x > box[0] && s.min_p.x < box[2] && + s.min_p.y > box[1] && s.min_p.y < box[3] && + s.max_p.x > box[0] && s.max_p.y < box[2] && + s.max_p.y > box[1] && s.max_p.y < box[3]) + score += boxes(i); + } + sum += score; + prior.insert({s.id, score}); + } + +// float max = boxes(0); +// float min = boxes(boxes.size() - 1); +// for (int i = 0; i < boxes.size(); i++) { +// float score = (boxes(i) - min) / (max - min); // normalized score +// // std::cout << "score: " << score << std::endl; +// for (int x = boxes[i][0] - 1; x < boxes[i][2]; x++) { +// for (int y = boxes[i][1] - 1; y < boxes[i][3]; y++) { +// likelihood.at(y, x) += score; +// } +// } +// } +// +// float sum = 0.f; +// for (const auto & s: segments) { +// cv::Mat mask = cv::Mat::zeros(image.size(), CV_8UC1); +// cv::rectangle(mask, s.min_p, s.max_p, cv::Scalar(1), CV_FILLED); +// float score = cv::mean(likelihood, mask)[0]; +// sum += score; +// prior.insert({s.id, score}); +// } + + for (auto & p: prior) { + p.second /= sum; + } + +// cv::normalize(likelihood, likelihood, 0.f, 1.f, cv::NORM_MINMAX); +// cv::namedWindow("Likelihood", cv::WINDOW_AUTOSIZE); +// cv::namedWindow("Image", cv::WINDOW_AUTOSIZE); +// cv::imshow("Likelihood", likelihood); +// cv::imshow("Image", image); +// cv::waitKey(); + + return prior; + } + +protected: + Objectness obj; +}; diff --git a/random/src/random_stopping_criterion.h b/random/src/random_stopping_criterion.h new file mode 100644 index 0000000..e7bbeb3 --- /dev/null +++ b/random/src/random_stopping_criterion.h @@ -0,0 +1,24 @@ +#pragma once + +#include "stopping_criterion.h" + +class RandomStoppingCriterion : public StoppingCriterion { +public: + RandomStoppingCriterion() : threshold(0.95f) {//threshold(dis(gen)) { + } + + bool stop(const cv::Mat & image, cv::Rect roi) { + return dis(gen) > threshold; + } + +protected: + float threshold; + + static std::random_device rd; + static std::mt19937 gen; + static std::uniform_real_distribution dis; +}; + +std::random_device RandomStoppingCriterion::rd; +std::mt19937 RandomStoppingCriterion::gen(rd()); +std::uniform_real_distribution RandomStoppingCriterion::dis(0.f, 1.f); diff --git a/random/src/segment.h b/random/src/segment.h new file mode 100644 index 0000000..ac77aef --- /dev/null +++ b/random/src/segment.h @@ -0,0 +1,166 @@ +#pragma once + +#include + +enum { + COLOR_SIMILARITY = 0x1, + TEXTURE_SIMILARITY = 0x2, + SIZE_SIMILARITY = 0x4, + BBOX_SIMILARITY = 0x8, +}; + +#define HIST_SIZE 25 +#define HIST_STEP int(255 / HIST_SIZE) + +class Segment +{ +public: + Segment() {} + Segment(int id, int nchannels, cv::Size size) : +#ifdef DEBUG + mask(cv::Mat::zeros(size, CV_8UC1)), +#endif + size(0), + im_size(size.area()), + im_size_cv(size), + min_p(size.width, size.height), + max_p(0, 0), + id(id), + nchannels(nchannels) + { + color_hist.assign(nchannels * HIST_SIZE, 0.f); + texture_hist.assign(HIST_SIZE, 0.f); + history.insert(id); + } + + template + void addPoint(cv::Vec<_Tp, n> col, uint8_t text, cv::Point point) { + size++; + + for (int i = 0; i < n; i++) + color_hist[i*HIST_SIZE + col[i] / HIST_STEP]++; + + texture_hist[text / HIST_STEP]++; + +#ifdef DEBUG + mask.at(point) = 1; +#endif + + min_p = cv::Point(std::min(point.x, min_p.x), std::min(point.y, min_p.y)); + max_p = cv::Point(std::max(point.x, max_p.x), std::max(point.y, max_p.y)); + } + + inline void addNeighbour(int n) { + neighbours.insert(n); + } + + void normalizeHistogram() { + cv::normalize(color_hist, color_hist, 1, 0, cv::NORM_L1); + cv::normalize(texture_hist, texture_hist, 1, 0, cv::NORM_L1); + } + + static float computeSimilarity(const Segment & a, const Segment & b, uint8_t flags) + { + if (a.im_size != b.im_size) { + std::cerr << "Error! Segments are from different images?" << std::endl; + return 0.f; + } + if (a.size == 0 || b.size == 0) { + std::cout << "Warning! Invalid segment." << std::endl; + return 0.f; + } + if (a.color_hist.size() != b.color_hist.size() || a.texture_hist.size() != b.texture_hist.size()) { + std::cout << "Error! Segments histograms don't match!" << std::endl; + return 0.f; + } + + float similarity = 0.f; + + if (flags & COLOR_SIMILARITY) { + for (uint64_t i = 0; i < a.color_hist.size(); i++) + similarity += std::min(a.color_hist[i], b.color_hist[i]); + } + + if (flags & TEXTURE_SIMILARITY) { + for (uint64_t i = 0; i < a.texture_hist.size(); i++) + similarity += std::min(a.texture_hist[i], b.texture_hist[i]); + } + + if (flags & SIZE_SIMILARITY) { + similarity += 1.f - float(a.size + b.size) / a.im_size; + } + + if (flags & BBOX_SIMILARITY) { + cv::Point min_p(std::min(a.min_p.x, b.min_p.x), std::min(a.min_p.y, b.min_p.y)); + cv::Point max_p(std::max(a.max_p.x, b.max_p.x), std::max(a.max_p.y, b.max_p.y)); + int bbox_size = (max_p.x - min_p.x) * (max_p.y - min_p.y); + similarity += 1.f - float(bbox_size - a.size - b.size) / a.im_size; + } + + return similarity; + } + + static Segment merge(const Segment & a, const Segment & b) { + Segment s(-1, a.nchannels, a.im_size_cv); + s.size = a.size + b.size; +#ifdef DEBUG + cv::bitwise_or(a.mask, b.mask, s.mask); +#endif + + if (a.color_hist.size() != b.color_hist.size() || a.texture_hist.size() != b.texture_hist.size()) + return s; + + for (uint64_t i = 0; i < a.color_hist.size(); i++) { + s.color_hist[i] = (a.color_hist[i] * a.size + b.size * b.color_hist[i]) / s.size; + } + + for (uint64_t i = 0; i < a.texture_hist.size(); i++) + s.texture_hist[i] = (a.size * a.texture_hist[i] + b.size * b.texture_hist[i]) / s.size; + //s.texture_hist = (a.size * a.texture_hist + b.size * b.texture_hist) / s.size; + + s.min_p = cv::Point(std::min(a.min_p.x, b.min_p.x), std::min(a.min_p.y, b.min_p.y)); + s.max_p = cv::Point(std::max(a.max_p.x, b.max_p.x), std::max(a.max_p.y, b.max_p.y)); + + s.history.insert(a.history.begin(), a.history.end()); + s.history.insert(b.history.begin(), b.history.end()); + + s.neighbours.insert(a.neighbours.begin(), a.neighbours.end()); + s.neighbours.insert(b.neighbours.begin(), b.neighbours.end()); + + for (const auto & h: s.history) + s.neighbours.erase(h); + + return s; + } + + friend std::ostream& operator<<(std::ostream &out, Segment & s) { + out << ""; + } + + inline cv::Rect bbox() { + return cv::Rect(min_p, max_p); + } + + inline bool empty() const { + return size == 0; + } + + std::set neighbours; + std::set history; +#ifdef DEBUG + cv::Mat mask; +#endif + int size; + int im_size; + cv::Size im_size_cv; + std::vector color_hist; + std::vector texture_hist; + cv::Point min_p; + cv::Point max_p; + int id; + int nchannels; +}; + diff --git a/random/src/selection_prior.h b/random/src/selection_prior.h new file mode 100644 index 0000000..6918842 --- /dev/null +++ b/random/src/selection_prior.h @@ -0,0 +1,37 @@ +#pragma once + +#include +#include +#include "segment.h" + +class SelectionPriorMap : public std::unordered_map { +public: + SelectionPriorMap() : std::unordered_map() {} + + uint32_t poll() { + float rnd = dis(gen); + float sum = 0.f; + for (auto p: *this) { + sum += p.second; + if (sum >= rnd) { + return p.first; + } + } + + return -1; + } + +protected: + static std::random_device rd; + static std::mt19937 gen; + static std::uniform_real_distribution dis; +}; + +std::random_device SelectionPriorMap::rd; +std::mt19937 SelectionPriorMap::gen(rd()); +std::uniform_real_distribution SelectionPriorMap::dis(0.f, 1.f); + +class SelectionPrior { +public: + virtual SelectionPriorMap computeSelectionPrior(const cv::Mat & image, const std::vector & segments) = 0; +}; diff --git a/random/src/stopping_criterion.h b/random/src/stopping_criterion.h new file mode 100644 index 0000000..e670ed7 --- /dev/null +++ b/random/src/stopping_criterion.h @@ -0,0 +1,8 @@ +#pragma once + +#include + +class StoppingCriterion { +public: + virtual bool stop(const cv::Mat & image, cv::Rect roi) = 0; +}; diff --git a/random/src/uniform_prior.h b/random/src/uniform_prior.h new file mode 100644 index 0000000..4d945ff --- /dev/null +++ b/random/src/uniform_prior.h @@ -0,0 +1,25 @@ +#pragma once + +#include "selection_prior.h" + +class UniformPrior : public SelectionPrior { +public: + virtual SelectionPriorMap computeSelectionPrior(const cv::Mat & image, const std::vector & segments) { + SelectionPriorMap prior; + + int sum = 0; + for (const auto s: segments) { + if (s.empty()) + continue; + + sum++; + prior.insert({s.id, 1}); + } + + for (auto & p: prior) { + p.second /= sum; + } + + return prior; + } +}; diff --git a/sande/CMakeLists.txt b/sande/CMakeLists.txt new file mode 100644 index 0000000..8457655 --- /dev/null +++ b/sande/CMakeLists.txt @@ -0,0 +1,36 @@ +cmake_minimum_required(VERSION 2.6) +project(segmentationtree) +set(PROJECT_NAME segmentationtree) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +# Boost +set(Boost_USE_STATIC_LIBS ON) +set(Boost_USE_MULTITHREADED ON) +find_package(Boost COMPONENTS python REQUIRED) + +find_package(OpenCV REQUIRED) +find_package(PythonLibs REQUIRED) +find_package(NumPy REQUIRED) + +set(CMAKE_CXX_FLAGS "-std=c++11 -O3") + +include_directories( + ${PYTHON_INCLUDE_DIRS} + ${NUMPY_INCLUDE_DIRS} +) + +add_library(selective_search MODULE src/main.cpp src/conversion.cpp) +set_target_properties(selective_search PROPERTIES PREFIX "") +target_link_libraries(selective_search + ${Boost_LIBRARIES} + ${OpenCV_LIBRARIES} + ${PYTHON_LIBRARIES} +) + +add_executable(test src/main.cpp src/conversion.cpp) +target_link_libraries(test + ${Boost_LIBRARIES} + ${OpenCV_LIBRARIES} + ${PYTHON_LIBRARIES} +) diff --git a/sande/cmake/FindNumPy.cmake b/sande/cmake/FindNumPy.cmake new file mode 100644 index 0000000..eafed16 --- /dev/null +++ b/sande/cmake/FindNumPy.cmake @@ -0,0 +1,102 @@ +# - Find the NumPy libraries +# This module finds if NumPy is installed, and sets the following variables +# indicating where it is. +# +# TODO: Update to provide the libraries and paths for linking npymath lib. +# +# NUMPY_FOUND - was NumPy found +# NUMPY_VERSION - the version of NumPy found as a string +# NUMPY_VERSION_MAJOR - the major version number of NumPy +# NUMPY_VERSION_MINOR - the minor version number of NumPy +# NUMPY_VERSION_PATCH - the patch version number of NumPy +# NUMPY_VERSION_DECIMAL - e.g. version 1.6.1 is 10601 +# NUMPY_INCLUDE_DIRS - path to the NumPy include files + +#============================================================================ +# Copyright 2012 Continuum Analytics, Inc. +# +# MIT License +# +# Permission is hereby granted, free of charge, to any person obtaining +# a copy of this software and associated documentation files +# (the "Software"), to deal in the Software without restriction, including +# without limitation the rights to use, copy, modify, merge, publish, +# distribute, sublicense, and/or sell copies of the Software, and to permit +# persons to whom the Software is furnished to do so, subject to +# the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR +# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +# OTHER DEALINGS IN THE SOFTWARE. +# +#============================================================================ + +# Finding NumPy involves calling the Python interpreter +if(NumPy_FIND_REQUIRED) + find_package(PythonInterp REQUIRED) +else() + find_package(PythonInterp) +endif() + +if(NOT PYTHONINTERP_FOUND) + set(NUMPY_FOUND FALSE) + return() +endif() + +execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-c" + "import numpy as n; print(n.__version__); print(n.get_include());" + RESULT_VARIABLE _NUMPY_SEARCH_SUCCESS + OUTPUT_VARIABLE _NUMPY_VALUES_OUTPUT + ERROR_VARIABLE _NUMPY_ERROR_VALUE + OUTPUT_STRIP_TRAILING_WHITESPACE) + +if(NOT _NUMPY_SEARCH_SUCCESS MATCHES 0) + if(NumPy_FIND_REQUIRED) + message(FATAL_ERROR + "NumPy import failure:\n${_NUMPY_ERROR_VALUE}") + endif() + set(NUMPY_FOUND FALSE) + return() +endif() + +# Convert the process output into a list +string(REGEX REPLACE ";" "\\\\;" _NUMPY_VALUES ${_NUMPY_VALUES_OUTPUT}) +string(REGEX REPLACE "\n" ";" _NUMPY_VALUES ${_NUMPY_VALUES}) +# Just in case there is unexpected output from the Python command. +list(GET _NUMPY_VALUES -2 NUMPY_VERSION) +list(GET _NUMPY_VALUES -1 NUMPY_INCLUDE_DIRS) + +string(REGEX MATCH "^[0-9]+\\.[0-9]+\\.[0-9]+" _VER_CHECK "${NUMPY_VERSION}") +if("${_VER_CHECK}" STREQUAL "") + # The output from Python was unexpected. Raise an error always + # here, because we found NumPy, but it appears to be corrupted somehow. + message(FATAL_ERROR + "Requested version and include path from NumPy, got instead:\n${_NUMPY_VALUES_OUTPUT}\n") + return() +endif() + +# Make sure all directory separators are '/' +string(REGEX REPLACE "\\\\" "/" NUMPY_INCLUDE_DIRS ${NUMPY_INCLUDE_DIRS}) + +# Get the major and minor version numbers +string(REGEX REPLACE "\\." ";" _NUMPY_VERSION_LIST ${NUMPY_VERSION}) +list(GET _NUMPY_VERSION_LIST 0 NUMPY_VERSION_MAJOR) +list(GET _NUMPY_VERSION_LIST 1 NUMPY_VERSION_MINOR) +list(GET _NUMPY_VERSION_LIST 2 NUMPY_VERSION_PATCH) +string(REGEX MATCH "[0-9]*" NUMPY_VERSION_PATCH ${NUMPY_VERSION_PATCH}) +math(EXPR NUMPY_VERSION_DECIMAL + "(${NUMPY_VERSION_MAJOR} * 10000) + (${NUMPY_VERSION_MINOR} * 100) + ${NUMPY_VERSION_PATCH}") + +find_package_message(NUMPY + "Found NumPy: version \"${NUMPY_VERSION}\" ${NUMPY_INCLUDE_DIRS}" + "${NUMPY_INCLUDE_DIRS}${NUMPY_VERSION}") + +set(NUMPY_FOUND TRUE) + diff --git a/sande/src/CVBoostConverter.hpp b/sande/src/CVBoostConverter.hpp new file mode 100644 index 0000000..5bbf456 --- /dev/null +++ b/sande/src/CVBoostConverter.hpp @@ -0,0 +1,423 @@ +/* + * CVBoostConverter.hpp + * + * Created on: Mar 20, 2014 + * Author: Gregory Kramida + * Copyright: (c) 2014 Gregory Kramida + * License: MIT + */ + +#ifndef CVBOOSTCONVERTER_HPP_ +#define CVBOOSTCONVERTER_HPP_ + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include +#include +#include +#include +#include +#include + +using namespace cv; +namespace bcvt{ +static PyObject* opencv_error = 0; + +class PyAllowThreads { +public: + PyAllowThreads() : + _state(PyEval_SaveThread()) { + } + ~PyAllowThreads() { + PyEval_RestoreThread(_state); + } +private: + PyThreadState* _state; +}; + +class PyEnsureGIL { +public: + PyEnsureGIL() : + _state(PyGILState_Ensure()) { + } + ~PyEnsureGIL() { + PyGILState_Release(_state); + } +private: + PyGILState_STATE _state; +}; +#define ERRWRAP2(expr) \ +try \ +{ \ + PyAllowThreads allowThreads; \ + expr; \ +} \ +catch (const cv::Exception &e) \ +{ \ + PyErr_SetString(opencv_error, e.what()); \ + return 0; \ +} + +class NumpyAllocator: public MatAllocator { +public: + NumpyAllocator() { + stdAllocator = Mat::getStdAllocator(); + } + ~NumpyAllocator() { + } + + UMatData* allocate(PyObject* o, int dims, const int* sizes, int type, + size_t* step) const { + UMatData* u = new UMatData(this); + u->data = u->origdata = (uchar*) PyArray_DATA((PyArrayObject*) o); + npy_intp* _strides = PyArray_STRIDES((PyArrayObject*) o); + for (int i = 0; i < dims - 1; i++) + step[i] = (size_t) _strides[i]; + step[dims - 1] = CV_ELEM_SIZE(type); + u->size = sizes[0] * step[0]; + u->userdata = o; + return u; + } + + UMatData* allocate(int dims0, const int* sizes, int type, void* data, + size_t* step, int flags, UMatUsageFlags usageFlags) const { + if (data != 0) { + CV_Error(Error::StsAssert, "The data should normally be NULL!"); + // probably this is safe to do in such extreme case + return stdAllocator->allocate(dims0, sizes, type, data, step, flags, + usageFlags); + } + PyEnsureGIL gil; + + int depth = CV_MAT_DEPTH(type); + int cn = CV_MAT_CN(type); + const int f = (int) (sizeof(size_t) / 8); + int typenum = + depth == CV_8U ? NPY_UBYTE : + depth == CV_8S ? NPY_BYTE : + depth == CV_16U ? NPY_USHORT : + depth == CV_16S ? NPY_SHORT : + depth == CV_32S ? NPY_INT : + depth == CV_32F ? NPY_FLOAT : + depth == CV_64F ? + NPY_DOUBLE : f * NPY_ULONGLONG + (f ^ 1) * NPY_UINT; + int i, dims = dims0; + cv::AutoBuffer _sizes(dims + 1); + for (i = 0; i < dims; i++) + _sizes[i] = sizes[i]; + if (cn > 1) + _sizes[dims++] = cn; + PyObject* o = PyArray_SimpleNew(dims, _sizes, typenum); + if (!o) + CV_Error_(Error::StsError, + ("The numpy array of typenum=%d, ndims=%d can not be created", typenum, dims)); + return allocate(o, dims0, sizes, type, step); + } + + bool allocate(UMatData* u, int accessFlags, + UMatUsageFlags usageFlags) const { + return stdAllocator->allocate(u, accessFlags, usageFlags); + } + + void deallocate(UMatData* u) const { + if (u) { + PyEnsureGIL gil; + PyObject* o = (PyObject*) u->userdata; + Py_XDECREF(o); + delete u; + } + } + + const MatAllocator* stdAllocator; +}; + +NumpyAllocator g_numpyAllocator; + +PyObject* fromMatToNDArray(const Mat& m) { + if (!m.data) + Py_RETURN_NONE; + Mat temp, + *p = (Mat*) &m; + if (!p->u || p->allocator != &g_numpyAllocator) { + temp.allocator = &g_numpyAllocator; + ERRWRAP2(m.copyTo(temp)); + p = &temp; + } + PyObject* o = (PyObject*) p->u->userdata; + Py_INCREF(o); + return o; +} + +static int failmsg(const char *fmt, ...) { + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +Mat fromNDArrayToMat(PyObject* o) { + cv::Mat m; + bool allowND = true; + if (!PyArray_Check(o)) { + failmsg("argument is not a numpy array"); + if (!m.data) + m.allocator = &g_numpyAllocator; + } else { + PyArrayObject* oarr = (PyArrayObject*) o; + + bool needcopy = false, needcast = false; + int typenum = PyArray_TYPE(oarr), new_typenum = typenum; + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : + typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT ? CV_32S : + typenum == NPY_INT32 ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if (type < 0) { + if (typenum == NPY_INT64 || typenum == NPY_UINT64 + || type == NPY_LONG) { + needcopy = needcast = true; + new_typenum = NPY_INT; + type = CV_32S; + } else { + failmsg("Argument data type is not supported"); + m.allocator = &g_numpyAllocator; + return m; + } + } + +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + + int ndims = PyArray_NDIM(oarr); + if (ndims >= CV_MAX_DIM) { + failmsg("Dimensionality of argument is too high"); + if (!m.data) + m.allocator = &g_numpyAllocator; + return m; + } + + int size[CV_MAX_DIM + 1]; + size_t step[CV_MAX_DIM + 1]; + size_t elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(oarr); + const npy_intp* _strides = PyArray_STRIDES(oarr); + bool ismultichannel = ndims == 3 && _sizes[2] <= CV_CN_MAX; + + for (int i = ndims - 1; i >= 0 && !needcopy; i--) { + // these checks handle cases of + // a) multi-dimensional (ndims > 2) arrays, as well as simpler 1- and 2-dimensional cases + // b) transposed arrays, where _strides[] elements go in non-descending order + // c) flipped arrays, where some of _strides[] elements are negative + if ((i == ndims - 1 && (size_t) _strides[i] != elemsize) + || (i < ndims - 1 && _strides[i] < _strides[i + 1])) + needcopy = true; + } + + if (ismultichannel && _strides[1] != (npy_intp) elemsize * _sizes[2]) + needcopy = true; + + if (needcopy) { + + if (needcast) { + o = PyArray_Cast(oarr, new_typenum); + oarr = (PyArrayObject*) o; + } else { + oarr = PyArray_GETCONTIGUOUS(oarr); + o = (PyObject*) oarr; + } + + _strides = PyArray_STRIDES(oarr); + } + + for (int i = 0; i < ndims; i++) { + size[i] = (int) _sizes[i]; + step[i] = (size_t) _strides[i]; + } + + // handle degenerate case + if (ndims == 0) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if (ismultichannel) { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + + if (ndims > 2 && !allowND) { + failmsg("%s has more than 2 dimensions"); + } else { + + m = Mat(ndims, size, type, PyArray_DATA(oarr), step); + m.u = g_numpyAllocator.allocate(o, ndims, size, type, step); + m.addref(); + + if (!needcopy) { + Py_INCREF(o); + } + } + m.allocator = &g_numpyAllocator; + } + return m; +} + +struct matToNDArrayBoostConverter { + static PyObject* convert(Mat const& m) { + if (!m.data) + Py_RETURN_NONE; + Mat *p = (Mat*) &m; + Mat temp; + if (!p->u || p->allocator != &g_numpyAllocator) { + temp.allocator = &g_numpyAllocator; + ERRWRAP2(m.copyTo(temp)); + p = &temp; + } + PyObject* o = (PyObject*) p->u->userdata; + return boost::python::incref(o); + } +}; + +struct matFromNDArrayBoostConverter { + + matFromNDArrayBoostConverter() { + boost::python::converter::registry::push_back(convertible, construct, + boost::python::type_id()); + } + + /// @brief Check if PyObject is an array and can be converted to OpenCV matrix. + static void* convertible(PyObject* object) { + if (!PyArray_Check(object)) { + return NULL; + } +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + PyArrayObject* oarr = (PyArrayObject*) object; + + int typenum = PyArray_TYPE(oarr); + if (typenum != NPY_INT64 && typenum != NPY_UINT64 && typenum != NPY_LONG + && typenum != NPY_UBYTE && typenum != NPY_BYTE + && typenum != NPY_USHORT && typenum != NPY_SHORT + && typenum != NPY_INT && typenum != NPY_INT32 + && typenum != NPY_FLOAT && typenum != NPY_DOUBLE) { + return NULL; + } + int ndims = PyArray_NDIM(oarr); //data type not supported + + if (ndims >= CV_MAX_DIM) { + return NULL; //too many dimensions + } + return object; + } + + /// @brief Construct a Mat from an NDArray object. + static void construct(PyObject* object, + boost::python::converter::rvalue_from_python_stage1_data* data) { + namespace python = boost::python; + // Object is a borrowed reference, so create a handle indicting it is + // borrowed for proper reference counting. + python::handle<> handle(python::borrowed(object)); + + // Obtain a handle to the memory block that the converter has allocated + // for the C++ type. + typedef python::converter::rvalue_from_python_storage storage_type; + void* storage = reinterpret_cast(data)->storage.bytes; + + // Allocate the C++ type into the converter's memory block, and assign + // its handle to the converter's convertible variable. The C++ + // container is populated by passing the begin and end iterators of + // the python object to the container's constructor. + PyArrayObject* oarr = (PyArrayObject*) object; + + bool needcopy = false, needcast = false; + int typenum = PyArray_TYPE(oarr), new_typenum = typenum; + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : + typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT ? CV_32S : + typenum == NPY_INT32 ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if (type < 0) { + needcopy = needcast = true; + new_typenum = NPY_INT; + type = CV_32S; + } + +#ifndef CV_MAX_DIM + const int CV_MAX_DIM = 32; +#endif + int ndims = PyArray_NDIM(oarr); + + int size[CV_MAX_DIM + 1]; + size_t step[CV_MAX_DIM + 1]; + size_t elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(oarr); + const npy_intp* _strides = PyArray_STRIDES(oarr); + bool ismultichannel = ndims == 3 && _sizes[2] <= CV_CN_MAX; + + for (int i = ndims - 1; i >= 0 && !needcopy; i--) { + // these checks handle cases of + // a) multi-dimensional (ndims > 2) arrays, as well as simpler 1- and 2-dimensional cases + // b) transposed arrays, where _strides[] elements go in non-descending order + // c) flipped arrays, where some of _strides[] elements are negative + if ((i == ndims - 1 && (size_t) _strides[i] != elemsize) + || (i < ndims - 1 && _strides[i] < _strides[i + 1])) + needcopy = true; + } + + if (ismultichannel && _strides[1] != (npy_intp) elemsize * _sizes[2]) + needcopy = true; + + if (needcopy) { + + if (needcast) { + object = PyArray_Cast(oarr, new_typenum); + oarr = (PyArrayObject*) object; + } else { + oarr = PyArray_GETCONTIGUOUS(oarr); + object = (PyObject*) oarr; + } + + _strides = PyArray_STRIDES(oarr); + } + + for (int i = 0; i < ndims; i++) { + size[i] = (int) _sizes[i]; + step[i] = (size_t) _strides[i]; + } + + // handle degenerate case + if (ndims == 0) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if (ismultichannel) { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + if (!needcopy) { + Py_INCREF(object); + } + cv::Mat* m = new (storage) cv::Mat(ndims, size, type, PyArray_DATA(oarr), step); + m->u = g_numpyAllocator.allocate(object, ndims, size, type, step); + m->allocator = &g_numpyAllocator; + m->addref(); + data->convertible = storage; + } +}; +} // end namespace bcvt +#endif /* CVBOOSTCONVERTER_HPP_ */ diff --git a/sande/src/adjacency.h b/sande/src/adjacency.h new file mode 100644 index 0000000..f99b2b9 --- /dev/null +++ b/sande/src/adjacency.h @@ -0,0 +1,29 @@ +#ifndef ADJACENCY_H_ +#define ADJACENCY_H_ + +#include + +class AdjacencyMatrix : public std::vector { +public: + AdjacencyMatrix(int n) { + this->n = n; + assign((n + 1)*n/2, false); + } + + reference get(int i, int j) { + return at(index(i, j)); + } + +protected: + int n; + + int index(int i, int j) { + int tmp = i > j ? j : i; + j = i > j ? i : j; + i = tmp; + return (n * i) + j - ((i * (i+1)) / 2); + } +}; + +#endif + diff --git a/sande/src/connection.h b/sande/src/connection.h new file mode 100644 index 0000000..801321d --- /dev/null +++ b/sande/src/connection.h @@ -0,0 +1,80 @@ +#ifndef CONNECTION_H_ +#define CONNECTION_H_ + +#include "segment.h" + +class Connection +{ +public: + Connection(Segment a_, Segment b_, uint8_t flags) : + a(a_.id), b(b_.id), similarity(exp(Segment::computeSimilarity(a_, b_, flags))) + { + } + + void merge(std::vector & connections, std::vector & segments, uint8_t flags, float & similarity_sum) { + Segment & a_ = segments[a]; + Segment & b_ = segments[b]; + + Segment s(segments.size(), a_.nchannels, a_.im_size_cv); + s.size = a_.size + b_.size; + cv::bitwise_or(a_.mask, b_.mask, s.mask); + + if (a_.color_hist.size() != b_.color_hist.size() || a_.texture_hist.size() != b_.texture_hist.size()) + return; + + for (uint64_t i = 0; i < a_.color_hist.size(); i++) { + s.color_hist[i] = (a_.color_hist[i] * a_.size + b_.size * b_.color_hist[i]) / s.size; + } + + for (uint64_t i = 0; i < a_.texture_hist.size(); i++) + s.texture_hist[i] = (a_.size * a_.texture_hist[i] + b_.size * b_.texture_hist[i]) / s.size; + //s.texture_hist = (a_.size * a_.texture_hist + b_.size * b_.texture_hist) / s.size; + + s.min_p = cv::Point(min(a_.min_p.x, b_.min_p.x), min(a_.min_p.y, b_.min_p.y)); + s.max_p = cv::Point(max(a_.max_p.x, b_.max_p.x), max(a_.max_p.y, b_.max_p.y)); + + for (auto n: a_.neighbours) { + if (n != b) + s.neighbours.insert(n); + } + for (auto n: b_.neighbours) { + if (n != a) + s.neighbours.insert(n); + } + + for (auto it = connections.begin(); it != connections.end(); ) { + if (it->a == a || it->a == b || it->b == a || it->b == b) { + similarity_sum -= it->similarity; + it = connections.erase(it); + } + else + it++; + } + + for (auto n: s.neighbours) { + Segment & neighbour = segments[n]; + neighbour.neighbours.erase(a); + neighbour.neighbours.erase(b); + neighbour.neighbours.insert(s.id); + Connection c(s, neighbour, flags); + connections.push_back(c); + similarity_sum += c.similarity; + } + + segments.push_back(s); + } + + friend std::ostream& operator<<(std::ostream &out, Connection & c) { + return out << ""; + } + + bool operator<(const Connection & rhs) const { + return rhs.similarity < similarity; + } + + int a; + int b; + float similarity; +}; + +#endif diff --git a/sande/src/conversion.cpp b/sande/src/conversion.cpp new file mode 100755 index 0000000..475c670 --- /dev/null +++ b/sande/src/conversion.cpp @@ -0,0 +1,231 @@ +# include "conversion.h" +/* + * The following conversion functions are taken/adapted from OpenCV's cv2.cpp file + * inside modules/python/src2 folder. + */ + +static void init() +{ + import_array(); +} + +static int failmsg(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class PyAllowThreads +{ +public: + PyAllowThreads() : _state(PyEval_SaveThread()) {} + ~PyAllowThreads() + { + PyEval_RestoreThread(_state); + } +private: + PyThreadState* _state; +}; + +class PyEnsureGIL +{ +public: + PyEnsureGIL() : _state(PyGILState_Ensure()) {} + ~PyEnsureGIL() + { + PyGILState_Release(_state); + } +private: + PyGILState_STATE _state; +}; + +using namespace cv; + +static PyObject* failmsgp(const char *fmt, ...) +{ + char str[1000]; + + va_list ap; + va_start(ap, fmt); + vsnprintf(str, sizeof(str), fmt, ap); + va_end(ap); + + PyErr_SetString(PyExc_TypeError, str); + return 0; +} + +class NumpyAllocator : public MatAllocator +{ +public: + NumpyAllocator() {} + ~NumpyAllocator() {} + + void allocate(int dims, const int* sizes, int type, int*& refcount, + uchar*& datastart, uchar*& data, size_t* step) + { + PyEnsureGIL gil; + + int depth = CV_MAT_DEPTH(type); + int cn = CV_MAT_CN(type); + const int f = (int)(sizeof(size_t)/8); + int typenum = depth == CV_8U ? NPY_UBYTE : depth == CV_8S ? NPY_BYTE : + depth == CV_16U ? NPY_USHORT : depth == CV_16S ? NPY_SHORT : + depth == CV_32S ? NPY_INT : depth == CV_32F ? NPY_FLOAT : + depth == CV_64F ? NPY_DOUBLE : f*NPY_ULONGLONG + (f^1)*NPY_UINT; + int i; + npy_intp _sizes[CV_MAX_DIM+1]; + for( i = 0; i < dims; i++ ) + { + _sizes[i] = sizes[i]; + } + + if( cn > 1 ) + { + _sizes[dims++] = cn; + } + + PyObject* o = PyArray_SimpleNew(dims, _sizes, typenum); + + if(!o) + { + CV_Error_(CV_StsError, ("The numpy array of typenum=%d, ndims=%d can not be created", typenum, dims)); + } + refcount = refcountFromPyObject(o); + + npy_intp* _strides = PyArray_STRIDES(o); + for( i = 0; i < dims - (cn > 1); i++ ) + step[i] = (size_t)_strides[i]; + datastart = data = (uchar*)PyArray_DATA(o); + } + + void deallocate(int* refcount, uchar*, uchar*) + { + PyEnsureGIL gil; + if( !refcount ) + return; + PyObject* o = pyObjectFromRefcount(refcount); + Py_INCREF(o); + Py_DECREF(o); + } +}; + +NumpyAllocator g_numpyAllocator; + +NDArrayConverter::NDArrayConverter() { init(); } + +void NDArrayConverter::init() +{ + import_array(); +} + +cv::Mat NDArrayConverter::toMat(const PyObject *o) +{ + cv::Mat m; + + if(!o || o == Py_None) + { + if( !m.data ) + m.allocator = &g_numpyAllocator; + } + + if( !PyArray_Check(o) ) + { + failmsg("toMat: Object is not a numpy array"); + } + + int typenum = PyArray_TYPE(o); + int type = typenum == NPY_UBYTE ? CV_8U : typenum == NPY_BYTE ? CV_8S : + typenum == NPY_USHORT ? CV_16U : typenum == NPY_SHORT ? CV_16S : + typenum == NPY_INT || typenum == NPY_LONG ? CV_32S : + typenum == NPY_FLOAT ? CV_32F : + typenum == NPY_DOUBLE ? CV_64F : -1; + + if( type < 0 ) + { + failmsg("toMat: Data type = %d is not supported", typenum); + } + + int ndims = PyArray_NDIM(o); + + if(ndims >= CV_MAX_DIM) + { + failmsg("toMat: Dimensionality (=%d) is too high", ndims); + } + + int size[CV_MAX_DIM+1]; + size_t step[CV_MAX_DIM+1], elemsize = CV_ELEM_SIZE1(type); + const npy_intp* _sizes = PyArray_DIMS(o); + const npy_intp* _strides = PyArray_STRIDES(o); + bool transposed = false; + + for(int i = 0; i < ndims; i++) + { + size[i] = (int)_sizes[i]; + step[i] = (size_t)_strides[i]; + } + + if( ndims == 0 || step[ndims-1] > elemsize ) { + size[ndims] = 1; + step[ndims] = elemsize; + ndims++; + } + + if( ndims >= 2 && step[0] < step[1] ) + { + std::swap(size[0], size[1]); + std::swap(step[0], step[1]); + transposed = true; + } + + if( ndims == 3 && size[2] <= CV_CN_MAX && step[1] == elemsize*size[2] ) + { + ndims--; + type |= CV_MAKETYPE(0, size[2]); + } + + if( ndims > 2) + { + failmsg("toMat: Object has more than 2 dimensions"); + } + + m = Mat(ndims, size, type, PyArray_DATA(o), step); + + if( m.data ) + { + m.refcount = refcountFromPyObject(o); + m.addref(); // protect the original numpy array from deallocation + // (since Mat destructor will decrement the reference counter) + }; + m.allocator = &g_numpyAllocator; + + if( transposed ) + { + Mat tmp; + tmp.allocator = &g_numpyAllocator; + transpose(m, tmp); + m = tmp; + } + return m; +} + +PyObject* NDArrayConverter::toNDArray(const cv::Mat& m) +{ + if( !m.data ) + Py_RETURN_NONE; + Mat temp, *p = (Mat*)&m; + if(!p->refcount || p->allocator != &g_numpyAllocator) + { + temp.allocator = &g_numpyAllocator; + m.copyTo(temp); + p = &temp; + } + p->addref(); + return pyObjectFromRefcount(p->refcount); +} diff --git a/sande/src/conversion.h b/sande/src/conversion.h new file mode 100755 index 0000000..3975e17 --- /dev/null +++ b/sande/src/conversion.h @@ -0,0 +1,60 @@ +# ifndef __COVERSION_OPENCV_H__ +# define __COVERSION_OPENCV_H__ + +#include +#include +#include +#include +#include "numpy/ndarrayobject.h" + +static PyObject* opencv_error = 0; + +static int failmsg(const char *fmt, ...); + +class PyAllowThreads; + +class PyEnsureGIL; + +#define ERRWRAP2(expr) \ +try \ +{ \ + PyAllowThreads allowThreads; \ + expr; \ +} \ +catch (const cv::Exception &e) \ +{ \ + PyErr_SetString(opencv_error, e.what()); \ + return 0; \ +} + +static PyObject* failmsgp(const char *fmt, ...); + +static size_t REFCOUNT_OFFSET = (size_t)&(((PyObject*)0)->ob_refcnt) + + (0x12345678 != *(const size_t*)"\x78\x56\x34\x12\0\0\0\0\0")*sizeof(int); + +static inline PyObject* pyObjectFromRefcount(const int* refcount) +{ + return (PyObject*)((size_t)refcount - REFCOUNT_OFFSET); +} + +static inline int* refcountFromPyObject(const PyObject* obj) +{ + return (int*)((size_t)obj + REFCOUNT_OFFSET); +} + + +class NumpyAllocator; + +enum { ARG_NONE = 0, ARG_MAT = 1, ARG_SCALAR = 2 }; + +class NDArrayConverter +{ +private: + void init(); +public: + NDArrayConverter(); + cv::Mat toMat(const PyObject* o); + PyObject* toNDArray(const cv::Mat& mat); +}; + +# endif diff --git a/sande/src/main.cpp b/sande/src/main.cpp new file mode 100644 index 0000000..12c5627 --- /dev/null +++ b/sande/src/main.cpp @@ -0,0 +1,192 @@ +//#include "CVBoostConverter.hpp" +#include +#include "conversion.h" +#include + +#include "connection.h" +#include "adjacency.h" + +using namespace boost::python; + +std::random_device rd_; +std::mt19937 gen_(rd_()); + +//#define DEBUG 1 + +cv::Mat get_bboxes_(cv::Mat image, cv::Mat seg, cv::Mat edge, uint8_t flags, uint8_t iterations) { + double max_id_; + cv::minMaxIdx(seg, nullptr, &max_id_); + int max_id = max_id_; + + std::vector segments; + std::vector connections; + segments.reserve(max_id); + int nchannels = image.channels(); + cv::Size size = image.size(); + for (int i = 0; i <= max_id; i++) { + segments.push_back(Segment(i, nchannels, size)); + } + + { + AdjacencyMatrix adjacency(max_id + 1); + for (int i = 0; i < image.rows; i++) { + for (int j = 0; j < image.cols; j++) { + cv::Point p(j, i); + uint16_t id = seg.at(p); + if (nchannels == 1) + segments[id].addPoint(image.at>(p), edge.at(p), p); + else if (nchannels == 3) + segments[id].addPoint(image.at(p), edge.at(p), p); + + if (i < image.rows - 1) { + uint16_t n = seg.at(i+1, j); + if (n != id && adjacency.get(id, n) == false) { + adjacency.get(id, n) = true; + segments[id].addNeighbour(n); + segments[n].addNeighbour(id); + } + } + + if (j < image.cols - 1) { + uint16_t n = seg.at(i, j+1); + if (n != id && adjacency.get(id, n) == false) { + adjacency.get(id, n) = true; + segments[id].addNeighbour(n); + segments[n].addNeighbour(id); + } + } + } + } + } + + float similarity_sum = 0.f; + { + AdjacencyMatrix adjacency(max_id + 1); + for (auto & s: segments) { + s.normalizeHistogram(); + for (auto n: s.neighbours) { + if (adjacency.get(s.id, n) == false) { + adjacency.get(s.id, n) = true; + Connection c(s, segments[n], flags); + connections.push_back(c); + similarity_sum += c.similarity; + } + } + } + } + + std::vector final_segments(segments); + bool stochastic = iterations > 1; + for (uint8_t i = 0; i < iterations; i++) { + std::vector connections_(connections); + std::vector segments_(segments); + float similarity_sum_ = similarity_sum; + while (connections_.size() != 0) { + std::uniform_real_distribution dis(0.f, similarity_sum_); + std::sort(connections_.begin(), connections_.end()); + float sum = 0.f; + float rnd = dis(gen_); + for (auto it = connections_.begin(); it != connections_.end(); it++) { + sum += it->similarity; + if (sum >= rnd || stochastic == false) { + similarity_sum_ -= it->similarity; + Connection c = *it; + connections_.erase(it); + c.merge(connections_, segments_, flags, similarity_sum_); + break; + } + } + } + + final_segments.insert(final_segments.end(), segments_.begin() + segments.size(), segments_.end()); + } + + cv::Mat bboxes; + for (auto s: final_segments) { + if (s.size == 0) + continue; + + cv::Mat bbox = cv::Mat(1, 4, CV_32SC1); + bbox.at(0) = s.min_p.x; + bbox.at(1) = s.min_p.y; + bbox.at(2) = s.max_p.x; + bbox.at(3) = s.max_p.y; + if (bboxes.empty()) + bboxes = bbox; + else + cv::vconcat(bboxes, bbox, bboxes); + } + +#ifdef DEBUG + cv::Mat gray; + cv::cvtColor(image, gray, cv::COLOR_BGR2GRAY); + cv::namedWindow("Segment", cv::WINDOW_NORMAL); + for (auto s: final_segments) { + cv::Mat display; + cv::addWeighted(s.mask * 255, 0.5, gray, 0.5, 0.0, display); + cv::rectangle(display, s.min_p, s.max_p, cv::Scalar(255)); + cv::imshow("Segment", display); + cv::waitKey(); + } +#endif + + return bboxes; +} + +PyObject * get_bboxes(PyObject * image_, PyObject * seg_, PyObject * edge_, uint8_t flags, uint8_t iterations) { + NDArrayConverter cvt; + cv::Mat image, seg, edge; + image = cvt.toMat(image_); + seg = cvt.toMat(seg_); + edge = cvt.toMat(edge_); + return cvt.toNDArray(get_bboxes_(image, seg, edge, flags, iterations)); +} + +static void init_ar() { + Py_Initialize(); + import_array(); +} + +BOOST_PYTHON_MODULE(selective_search) { + init_ar(); + + //initialize converters + /*to_python_converter(); + bcvt::matFromNDArrayBoostConverter();*/ + + def("get_bboxes", get_bboxes); +} + +int main(int argc, char * argv[]) { + if (argc != 5) { + std::cout << "Usage: " << argv[0] << " " << std::endl; + return 0; + } + + cv::Mat image = cv::imread(argv[1]); + cv::Mat seg = cv::imread(argv[2], cv::IMREAD_UNCHANGED); + cv::Mat edge = cv::imread(argv[3], cv::IMREAD_UNCHANGED); + uint8_t iterations = atoi(argv[4]); + +// cv::namedWindow("Image", cv::WINDOW_NORMAL); +// cv::imshow("Image", (seg == 338) * 255); +// cv::waitKey(); + + std::clock_t begin = std::clock(); + cv::Mat bboxes = get_bboxes_(image, seg, edge, COLOR_SIMILARITY | TEXTURE_SIMILARITY | SIZE_SIMILARITY | BBOX_SIMILARITY, iterations); + std::clock_t end = std::clock(); + double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC; + std::cout << "Times passed in seconds: " << elapsed_secs << std::endl; + +// for (int i = 0; i < bboxes.rows; i++) { +// cv::Point p1(bboxes.at(i,0), bboxes.at(i,1)); +// cv::Point p2(bboxes.at(i,2), bboxes.at(i,3)); +// cv::Mat tmp; +// image.copyTo(tmp); +// cv::rectangle(tmp, p1, p2, cv::Scalar(255, 255, 0)); +// cv::imshow("Image", tmp); +// cv::waitKey(); +// } + return 0; +} diff --git a/sande/src/segment.h b/sande/src/segment.h new file mode 100644 index 0000000..7f99404 --- /dev/null +++ b/sande/src/segment.h @@ -0,0 +1,127 @@ +#ifndef SEGMENT_H_ +#define SEGMENT_H_ + +#include + +enum { + COLOR_SIMILARITY = 0x1, + TEXTURE_SIMILARITY = 0x2, + SIZE_SIMILARITY = 0x4, + BBOX_SIMILARITY = 0x8, +}; + +const cv::Mat se = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3,3)); + +#define min(a, b) (a < b ? a : b) +#define max(a, b) (a > b ? a : b) + +#define HIST_SIZE 25 +#define HIST_STEP int(255 / HIST_SIZE) + +class Segment +{ +public: + Segment(int id, int nchannels, cv::Size size) : + mask(cv::Mat::zeros(size, CV_8UC1)), + size(0), + im_size(size.area()), + im_size_cv(size), + min_p(size.width, size.height), + max_p(0, 0), + id(id), + nchannels(nchannels) + { + color_hist.assign(nchannels * HIST_SIZE, 0.f); + texture_hist.assign(HIST_SIZE, 0.f); + } + + template + void addPoint(cv::Vec<_Tp, n> col, uint8_t text, cv::Point point) { + size++; + + for (int i = 0; i < n; i++) + color_hist[i*HIST_SIZE + col[i] / HIST_STEP]++; + + texture_hist[text / HIST_STEP]++; + + mask.at(point) = 1; + + min_p = cv::Point(min(point.x, min_p.x), min(point.y, min_p.y)); + max_p = cv::Point(max(point.x, max_p.x), max(point.y, max_p.y)); + } + + inline void addNeighbour(int n) { + neighbours.insert(n); + } + + void normalizeHistogram() { + cv::normalize(color_hist, color_hist, 1, 0, cv::NORM_L1); + cv::normalize(texture_hist, texture_hist, 1, 0, cv::NORM_L1); + } + + static float computeSimilarity(Segment a, Segment b, uint8_t flags) + { + if (a.im_size != b.im_size) { + std::cerr << "Error! Segments are from different images?" << std::endl; + return 0.f; + } + if (a.size == 0 || b.size == 0) { + std::cout << "Warning! Invalid segment." << std::endl; + return 0.f; + } + if (a.color_hist.size() != b.color_hist.size() || a.texture_hist.size() != b.texture_hist.size()) { + std::cout << "Error! Segments histograms don't match!" << std::endl; + return 0.f; + } + + float similarity = 0.f; + + if (flags & COLOR_SIMILARITY) { + for (uint64_t i = 0; i < a.color_hist.size(); i++) + similarity += min(a.color_hist[i], b.color_hist[i]); + } + + if (flags & TEXTURE_SIMILARITY) { + for (uint64_t i = 0; i < a.texture_hist.size(); i++) + similarity += min(a.texture_hist[i], b.texture_hist[i]); + } + + if (flags & SIZE_SIMILARITY) { + similarity += 1.f - float(a.size + b.size) / a.im_size; + } + + if (flags & BBOX_SIMILARITY) { + cv::Point min_p(min(a.min_p.x, b.min_p.x), min(a.min_p.y, b.min_p.y)); + cv::Point max_p(max(a.max_p.x, b.max_p.x), max(a.max_p.y, b.max_p.y)); + int bbox_size = (max_p.x - min_p.x) * (max_p.y - min_p.y); + similarity += 1.f - float(bbox_size - a.size - b.size) / a.im_size; + } + + return similarity; + } + + friend std::ostream& operator<<(std::ostream &out, Segment & s) { + out << ""; + } + + inline cv::Rect bbox() { + return cv::Rect(min_p, max_p); + } + + std::set neighbours; + cv::Mat mask; + int size; + int im_size; + cv::Size im_size_cv; + std::vector color_hist; + std::vector texture_hist; + cv::Point min_p; + cv::Point max_p; + int id; + int nchannels; +}; + +#endif diff --git a/tree/src/adjacency.h b/tree/src/adjacency.h new file mode 100644 index 0000000..865b937 --- /dev/null +++ b/tree/src/adjacency.h @@ -0,0 +1,28 @@ +#ifndef ADJACENCY_H_ +#define ADJACENCY_H_ + +class AdjacencyMatrix : public std::vector { +public: + AdjacencyMatrix(int n) { + this->n = n; + resize((n + 1)*n/2); + } + + Edge& get(int i, int j) { + int tmp = i > j ? j : i; + j = i > j ? i : j; + i = tmp; + int ind = index(i, j); + return operator[](ind); + } + +protected: + int n; + + int index(int i, int j) { + return (n * i) + j - ((i * (i+1)) / 2); + } +}; + +#endif + diff --git a/tree/src/tree.cpp b/tree/src/tree.cpp index 56f5b7c..bc68b18 100644 --- a/tree/src/tree.cpp +++ b/tree/src/tree.cpp @@ -1,37 +1,14 @@ -#define DEBUG 1 +//#define DEBUG #include "CVBoostConverter.hpp" #include #include #include #include "edge.h" +#include "adjacency.h" using namespace boost::python; - -class AdjacencyMatrix : public std::vector { -public: - AdjacencyMatrix(int n) { - this->n = n; - resize((n + 1)*n/2); - } - - Edge& get(int i, int j) { - int tmp = i > j ? j : i; - j = i > j ? i : j; - i = tmp; - int ind = index(i, j); - return operator[](ind); - } - -protected: - int n; - - int index(int i, int j) { - return (n * i) + j - ((i * (i+1)) / 2); - } -}; - void updateSurface(cv::Point p, int s, std::unordered_map> & surfaces, cv::Size size) { cv::Point p1, p2; auto it = surfaces.find(s); @@ -51,12 +28,13 @@ void updateSurface(cv::Point p, int s, std::unordered_map> surfaces; @@ -65,12 +43,16 @@ cv::Mat makeTree(cv::Mat s, cv::Mat e) { for (int i = 0; i < nonzero.total(); i++) { cv::Point p = nonzero.at(i); - if (s.at(p.y-1, p.x) > 0 && s.at(p.y+1, p.x) > 0 && - s.at(p.y-1, p.x) != s.at(p.y+1, p.x)) { - edges.at(p.y, p.x) = 255; + // skip edge cases + if (p.x == 0 || p.x == s.cols - 1 || + p.y == 0 || p.y == s.rows - 1) + continue; - int a = s.at(p.y-1, p.x); - int b = s.at(p.y+1, p.x); + // up/down are different clusters? + int a = s.at(p.y-1, p.x); + int b = s.at(p.y+1, p.x); + if (a > 0 && b > 0 && a != b) { + //edges.at(p.y, p.x) = 255; updateSurface(p, a, surfaces, s.size()); updateSurface(p, b, surfaces, s.size()); @@ -79,19 +61,20 @@ cv::Mat makeTree(cv::Mat s, cv::Mat e) { c.edges.push_back(e.at(p)); c.updateBox(p); } - else if (s.at(p.y, p.x - 1) > 0 && s.at(p.y, p.x + 1) > 0 && - s.at(p.y, p.x - 1) != s.at(p.y, p.x + 1)) { - edges.at(p.y, p.x) = 255; - - int a = s.at(p.y, p.x-1); - int b = s.at(p.y, p.x+1); - - updateSurface(p, a, surfaces, s.size()); - updateSurface(p, b, surfaces, s.size()); - - Edge & c = adjacency.get(a, b); - c.edges.push_back(e.at(p)); - c.updateBox(p); + // left/right are different clusters? + else { + a = s.at(p.y, p.x - 1); + b = s.at(p.y, p.x + 1); + if (a > 0 && b > 0 && a != b) { + //edges.at(p.y, p.x) = 255; + + updateSurface(p, a, surfaces, s.size()); + updateSurface(p, b, surfaces, s.size()); + + Edge & c = adjacency.get(a, b); + c.edges.push_back(e.at(p)); + c.updateBox(p); + } } } @@ -214,22 +197,20 @@ cv::Mat makeTree(cv::Mat s, cv::Mat e) { return result; } -static void init_ar() -{ - Py_Initialize(); - import_array(); +static void init_ar() { + Py_Initialize(); + import_array(); } -BOOST_PYTHON_MODULE(segmentation_tree) -{ - init_ar(); +BOOST_PYTHON_MODULE(segmentation_tree) { + init_ar(); //initialize converters to_python_converter(); bcvt::matFromNDArrayBoostConverter(); - def("tree", makeTree); + def("tree", makeTree); } int main(int argc, char * argv[]) { @@ -254,8 +235,11 @@ int main(int argc, char * argv[]) { cv::imshow("Segmentation", segmentation); cv::imshow("Edge", edge); + cv::Mat tree; std::clock_t begin = std::clock(); - cv::Mat tree = makeTree(segmentation, edge); + for (int i = 0; i < 1; i++) { + tree = makeTree(segmentation, edge); + } std::clock_t end = std::clock(); std::cout << tree.rows << " surfaces detected." << std::endl; double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;