generic-library/opencv
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

removed optim module; moved its functionality to core and photo modules; moved drawing functions from core to imgproc. Removed FilterEngine etc. from public API

Browse Source
This commit is contained in:
Vadim Pisarevsky
2014-08-14 12:50:07 +04:00
parent a602185fb6
commit 257463719b
44 changed files with 5886 additions and 3186 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
modules/core/doc/core.rst
View File

@@ -17,3 +17,4 @@ core. The Core Functionality
utility_and_system_functions_and_macros
opengl_interop
ipp_async_converters
optim

634
modules/core/doc/drawing_functions.rst
View File

@@ -1,634 +0,0 @@
Drawing Functions
=================
.. highlight:: cpp
Drawing functions work with matrices/images of arbitrary depth.
The boundaries of the shapes can be rendered with antialiasing (implemented only for 8-bit images for now).
All the functions include the parameter ``color`` that uses an RGB value (that may be constructed
with the ``Scalar`` constructor
) for color
images and brightness for grayscale images. For color images, the channel ordering
is normally *Blue, Green, Red*.
This is what :ocv:func:`imshow`, :ocv:func:`imread`, and :ocv:func:`imwrite` expect.
So, if you form a color using the
``Scalar`` constructor, it should look like:
.. math::
\texttt{Scalar} (blue \_ component, green \_ component, red \_ component[, alpha \_ component])
If you are using your own image rendering and I/O functions, you can use any channel ordering. The drawing functions process each channel independently and do not depend on the channel order or even on the used color space. The whole image can be converted from BGR to RGB or to a different color space using
:ocv:func:`cvtColor` .
If a drawn figure is partially or completely outside the image, the drawing functions clip it. Also, many drawing functions can handle pixel coordinates specified with sub-pixel accuracy. This means that the coordinates can be passed as fixed-point numbers encoded as integers. The number of fractional bits is specified by the ``shift`` parameter and the real point coordinates are calculated as
:math:`\texttt{Point}(x,y)\rightarrow\texttt{Point2f}(x*2^{-shift},y*2^{-shift})` . This feature is especially effective when rendering antialiased shapes.
.. note:: The functions do not support alpha-transparency when the target image is 4-channel. In this case, the ``color[3]`` is simply copied to the repainted pixels. Thus, if you want to paint semi-transparent shapes, you can paint them in a separate buffer and then blend it with the main image.
.. note::
* An example on using variate drawing functions like line, rectangle, ... can be found at opencv_source_code/samples/cpp/drawing.cpp
circle
----------
Draws a circle.
.. ocv:function:: void circle( InputOutputArray img, Point center, int radius, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:pyfunction:: cv2.circle(img, center, radius, color[, thickness[, lineType[, shift]]]) -> img
.. ocv:cfunction:: void cvCircle( CvArr* img, CvPoint center, int radius, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
:param img: Image where the circle is drawn.
:param center: Center of the circle.
:param radius: Radius of the circle.
:param color: Circle color.
:param thickness: Thickness of the circle outline, if positive. Negative thickness means that a filled circle is to be drawn.
:param lineType: Type of the circle boundary. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the coordinates of the center and in the radius value.
The function ``circle`` draws a simple or filled circle with a given center and radius.
clipLine
------------
Clips the line against the image rectangle.
.. ocv:function:: bool clipLine(Size imgSize, Point& pt1, Point& pt2)
.. ocv:function:: bool clipLine(Rect imgRect, Point& pt1, Point& pt2)
.. ocv:pyfunction:: cv2.clipLine(imgRect, pt1, pt2) -> retval, pt1, pt2
.. ocv:cfunction:: int cvClipLine( CvSize img_size, CvPoint* pt1, CvPoint* pt2 )
:param imgSize: Image size. The image rectangle is ``Rect(0, 0, imgSize.width, imgSize.height)`` .
:param imgRect: Image rectangle.
:param pt1: First line point.
:param pt2: Second line point.
The functions ``clipLine`` calculate a part of the line segment that is entirely within the specified rectangle.
They return ``false`` if the line segment is completely outside the rectangle. Otherwise, they return ``true`` .
ellipse
-----------
Draws a simple or thick elliptic arc or fills an ellipse sector.
.. ocv:function:: void ellipse( InputOutputArray img, Point center, Size axes, double angle, double startAngle, double endAngle, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:function:: void ellipse( InputOutputArray img, const RotatedRect& box, const Scalar& color, int thickness=1, int lineType=LINE_8 )
.. ocv:pyfunction:: cv2.ellipse(img, center, axes, angle, startAngle, endAngle, color[, thickness[, lineType[, shift]]]) -> img
.. ocv:pyfunction:: cv2.ellipse(img, box, color[, thickness[, lineType]]) -> img
.. ocv:cfunction:: void cvEllipse( CvArr* img, CvPoint center, CvSize axes, double angle, double start_angle, double end_angle, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
.. ocv:cfunction:: void cvEllipseBox( CvArr* img, CvBox2D box, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
:param img: Image.
:param center: Center of the ellipse.
:param axes: Half of the size of the ellipse main axes.
:param angle: Ellipse rotation angle in degrees.
:param startAngle: Starting angle of the elliptic arc in degrees.
:param endAngle: Ending angle of the elliptic arc in degrees.
:param box: Alternative ellipse representation via :ocv:class:`RotatedRect` or ``CvBox2D``. This means that the function draws an ellipse inscribed in the rotated rectangle.
:param color: Ellipse color.
:param thickness: Thickness of the ellipse arc outline, if positive. Otherwise, this indicates that a filled ellipse sector is to be drawn.
:param lineType: Type of the ellipse boundary. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the coordinates of the center and values of axes.
The functions ``ellipse`` with less parameters draw an ellipse outline, a filled ellipse, an elliptic arc, or a filled ellipse sector.
A piecewise-linear curve is used to approximate the elliptic arc boundary. If you need more control of the ellipse rendering, you can retrieve the curve using
:ocv:func:`ellipse2Poly` and then render it with
:ocv:func:`polylines` or fill it with
:ocv:func:`fillPoly` . If you use the first variant of the function and want to draw the whole ellipse, not an arc, pass ``startAngle=0`` and ``endAngle=360`` . The figure below explains the meaning of the parameters.
**Figure 1. Parameters of Elliptic Arc**
.. image:: pics/ellipse.png
ellipse2Poly
----------------
Approximates an elliptic arc with a polyline.
.. ocv:function:: void ellipse2Poly( Point center, Size axes, int angle, int arcStart, int arcEnd, int delta, vector<Point>& pts )
.. ocv:pyfunction:: cv2.ellipse2Poly(center, axes, angle, arcStart, arcEnd, delta) -> pts
:param center: Center of the arc.
:param axes: Half of the size of the ellipse main axes. See the :ocv:func:`ellipse` for details.
:param angle: Rotation angle of the ellipse in degrees. See the :ocv:func:`ellipse` for details.
:param arcStart: Starting angle of the elliptic arc in degrees.
:param arcEnd: Ending angle of the elliptic arc in degrees.
:param delta: Angle between the subsequent polyline vertices. It defines the approximation accuracy.
:param pts: Output vector of polyline vertices.
The function ``ellipse2Poly`` computes the vertices of a polyline that approximates the specified elliptic arc. It is used by
:ocv:func:`ellipse` .
fillConvexPoly
------------------
Fills a convex polygon.
.. ocv:function:: void fillConvexPoly( Mat& img, const Point* pts, int npts, const Scalar& color, int lineType=LINE_8, int shift=0 )
.. ocv:function:: void fillConvexPoly( InputOutputArray img, InputArray points, const Scalar& color, int lineType=LINE_8, int shift=0 )
.. ocv:pyfunction:: cv2.fillConvexPoly(img, points, color[, lineType[, shift]]) -> img
.. ocv:cfunction:: void cvFillConvexPoly( CvArr* img, const CvPoint* pts, int npts, CvScalar color, int line_type=8, int shift=0 )
:param img: Image.
:param pts: Polygon vertices.
:param npts: Number of polygon vertices.
:param color: Polygon color.
:param lineType: Type of the polygon boundaries. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the vertex coordinates.
The function ``fillConvexPoly`` draws a filled convex polygon.
This function is much faster than the function ``fillPoly`` . It can fill not only convex polygons but any monotonic polygon without self-intersections,
that is, a polygon whose contour intersects every horizontal line (scan line) twice at the most (though, its top-most and/or the bottom edge could be horizontal).
fillPoly
------------
Fills the area bounded by one or more polygons.
.. ocv:function:: void fillPoly( Mat& img, const Point** pts, const int* npts, int ncontours, const Scalar& color, int lineType=LINE_8, int shift=0, Point offset=Point() )
.. ocv:function:: void fillPoly( InputOutputArray img, InputArrayOfArrays pts, const Scalar& color, int lineType=LINE_8, int shift=0, Point offset=Point() )
.. ocv:pyfunction:: cv2.fillPoly(img, pts, color[, lineType[, shift[, offset]]]) -> img
.. ocv:cfunction:: void cvFillPoly( CvArr* img, CvPoint** pts, const int* npts, int contours, CvScalar color, int line_type=8, int shift=0 )
:param img: Image.
:param pts: Array of polygons where each polygon is represented as an array of points.
:param npts: Array of polygon vertex counters.
:param ncontours: Number of contours that bind the filled region.
:param color: Polygon color.
:param lineType: Type of the polygon boundaries. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the vertex coordinates.
:param offset: Optional offset of all points of the contours.
The function ``fillPoly`` fills an area bounded by several polygonal contours. The function can fill complex areas, for example,
areas with holes, contours with self-intersections (some of their parts), and so forth.
getTextSize
---------------
Calculates the width and height of a text string.
.. ocv:function:: Size getTextSize(const String& text, int fontFace, double fontScale, int thickness, int* baseLine)
.. ocv:pyfunction:: cv2.getTextSize(text, fontFace, fontScale, thickness) -> retval, baseLine
.. ocv:cfunction:: void cvGetTextSize( const char* text_string, const CvFont* font, CvSize* text_size, int* baseline )
:param text: Input text string.
:param text_string: Input text string in C format.
:param fontFace: Font to use. See the :ocv:func:`putText` for details.
:param fontScale: Font scale. See the :ocv:func:`putText` for details.
:param thickness: Thickness of lines used to render the text. See :ocv:func:`putText` for details.
:param baseLine: Output parameter - y-coordinate of the baseline relative to the bottom-most text point.
:param baseline: Output parameter - y-coordinate of the baseline relative to the bottom-most text point.
:param font: Font description in terms of old C API.
:param text_size: Output parameter - The size of a box that contains the specified text.
The function ``getTextSize`` calculates and returns the size of a box that contains the specified text.
That is, the following code renders some text, the tight box surrounding it, and the baseline: ::
String text = "Funny text inside the box";
int fontFace = FONT_HERSHEY_SCRIPT_SIMPLEX;
double fontScale = 2;
int thickness = 3;
Mat img(600, 800, CV_8UC3, Scalar::all(0));
int baseline=0;
Size textSize = getTextSize(text, fontFace,
fontScale, thickness, &baseline);
baseline += thickness;
// center the text
Point textOrg((img.cols - textSize.width)/2,
(img.rows + textSize.height)/2);
// draw the box
rectangle(img, textOrg + Point(0, baseline),
textOrg + Point(textSize.width, -textSize.height),
Scalar(0,0,255));
// ... and the baseline first
line(img, textOrg + Point(0, thickness),
textOrg + Point(textSize.width, thickness),
Scalar(0, 0, 255));
// then put the text itself
putText(img, text, textOrg, fontFace, fontScale,
Scalar::all(255), thickness, 8);
InitFont
--------
Initializes font structure (OpenCV 1.x API).
.. ocv:cfunction:: void cvInitFont( CvFont* font, int font_face, double hscale, double vscale, double shear=0, int thickness=1, int line_type=8 )
:param font: Pointer to the font structure initialized by the function
:param font_face: Font name identifier. Only a subset of Hershey fonts http://sources.isc.org/utils/misc/hershey-font.txt are supported now:
* **CV_FONT_HERSHEY_SIMPLEX** normal size sans-serif font
* **CV_FONT_HERSHEY_PLAIN** small size sans-serif font
* **CV_FONT_HERSHEY_DUPLEX** normal size sans-serif font (more complex than ``CV_FONT_HERSHEY_SIMPLEX`` )
* **CV_FONT_HERSHEY_COMPLEX** normal size serif font
* **CV_FONT_HERSHEY_TRIPLEX** normal size serif font (more complex than ``CV_FONT_HERSHEY_COMPLEX`` )
* **CV_FONT_HERSHEY_COMPLEX_SMALL** smaller version of ``CV_FONT_HERSHEY_COMPLEX``
* **CV_FONT_HERSHEY_SCRIPT_SIMPLEX** hand-writing style font
* **CV_FONT_HERSHEY_SCRIPT_COMPLEX** more complex variant of ``CV_FONT_HERSHEY_SCRIPT_SIMPLEX``
The parameter can be composited from one of the values above and an optional ``CV_FONT_ITALIC`` flag, which indicates italic or oblique font.
:param hscale: Horizontal scale. If equal to ``1.0f`` , the characters have the original width depending on the font type. If equal to ``0.5f`` , the characters are of half the original width.
:param vscale: Vertical scale. If equal to ``1.0f`` , the characters have the original height depending on the font type. If equal to ``0.5f`` , the characters are of half the original height.
:param shear: Approximate tangent of the character slope relative to the vertical line. A zero value means a non-italic font, ``1.0f`` means about a 45 degree slope, etc.
:param thickness: Thickness of the text strokes
:param line_type: Type of the strokes, see :ocv:func:`line` description
The function initializes the font structure that can be passed to text rendering functions.
.. seealso:: :ocv:cfunc:`PutText`
.. _Line:
line
--------
Draws a line segment connecting two points.
.. ocv:function:: void line( InputOutputArray img, Point pt1, Point pt2, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:pyfunction:: cv2.line(img, pt1, pt2, color[, thickness[, lineType[, shift]]]) -> img
.. ocv:cfunction:: void cvLine( CvArr* img, CvPoint pt1, CvPoint pt2, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
:param img: Image.
:param pt1: First point of the line segment.
:param pt2: Second point of the line segment.
:param color: Line color.
:param thickness: Line thickness.
:param lineType: Type of the line:
* **LINE_8** (or omitted) - 8-connected line.
* **LINE_4** - 4-connected line.
* **LINE_AA** - antialiased line.
:param shift: Number of fractional bits in the point coordinates.
The function ``line`` draws the line segment between ``pt1`` and ``pt2`` points in the image. The line is clipped by the image boundaries. For non-antialiased lines with integer coordinates, the 8-connected or 4-connected Bresenham algorithm is used. Thick lines are drawn with rounding endings.
Antialiased lines are drawn using Gaussian filtering.
arrowedLine
----------------
Draws a arrow segment pointing from the first point to the second one.
.. ocv:function:: void arrowedLine(InputOutputArray img, Point pt1, Point pt2, const Scalar& color, int thickness=1, int lineType=8, int shift=0, double tipLength=0.1)
:param img: Image.
:param pt1: The point the arrow starts from.
:param pt2: The point the arrow points to.
:param color: Line color.
:param thickness: Line thickness.
:param lineType: Type of the line:
* **8** (or omitted) - 8-connected line.
* **4** - 4-connected line.
* **CV_AA** - antialiased line.
:param shift: Number of fractional bits in the point coordinates.
:param tipLength: The length of the arrow tip in relation to the arrow length
The function ``arrowedLine`` draws an arrow between ``pt1`` and ``pt2`` points in the image. See also :ocv:func:`line`.
LineIterator
------------
.. ocv:class:: LineIterator
Class for iterating pixels on a raster line. ::
class LineIterator
{
public:
// creates iterators for the line connecting pt1 and pt2
// the line will be clipped on the image boundaries
// the line is 8-connected or 4-connected
// If leftToRight=true, then the iteration is always done
// from the left-most point to the right most,
// not to depend on the ordering of pt1 and pt2 parameters
LineIterator(const Mat& img, Point pt1, Point pt2,
int connectivity=8, bool leftToRight=false);
// returns pointer to the current line pixel
uchar* operator *();
// move the iterator to the next pixel
LineIterator& operator ++();
LineIterator operator ++(int);
Point pos() const;
// internal state of the iterator
uchar* ptr;
int err, count;
int minusDelta, plusDelta;
int minusStep, plusStep;
};
The class ``LineIterator`` is used to get each pixel of a raster line. It can be treated as versatile implementation of the Bresenham algorithm where you can stop at each pixel and do some extra processing, for example, grab pixel values along the line or draw a line with an effect (for example, with XOR operation).
The number of pixels along the line is stored in ``LineIterator::count`` . The method ``LineIterator::pos`` returns the current position in the image ::
// grabs pixels along the line (pt1, pt2)
// from 8-bit 3-channel image to the buffer
LineIterator it(img, pt1, pt2, 8);
LineIterator it2 = it;
vector<Vec3b> buf(it.count);
for(int i = 0; i < it.count; i++, ++it)
buf[i] = *(const Vec3b)*it;
// alternative way of iterating through the line
for(int i = 0; i < it2.count; i++, ++it2)
{
Vec3b val = img.at<Vec3b>(it2.pos());
CV_Assert(buf[i] == val);
}
rectangle
-------------
Draws a simple, thick, or filled up-right rectangle.
.. ocv:function:: void rectangle( InputOutputArray img, Point pt1, Point pt2, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:function:: void rectangle( Mat& img, Rect rec, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:pyfunction:: cv2.rectangle(img, pt1, pt2, color[, thickness[, lineType[, shift]]]) -> img
.. ocv:cfunction:: void cvRectangle( CvArr* img, CvPoint pt1, CvPoint pt2, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
:param img: Image.
:param pt1: Vertex of the rectangle.
:param pt2: Vertex of the rectangle opposite to ``pt1`` .
:param rec: Alternative specification of the drawn rectangle.
:param color: Rectangle color or brightness (grayscale image).
:param thickness: Thickness of lines that make up the rectangle. Negative values, like ``CV_FILLED`` , mean that the function has to draw a filled rectangle.
:param lineType: Type of the line. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the point coordinates.
The function ``rectangle`` draws a rectangle outline or a filled rectangle whose two opposite corners are ``pt1`` and ``pt2``, or ``r.tl()`` and ``r.br()-Point(1,1)``.
polylines
-------------
Draws several polygonal curves.
.. ocv:function:: void polylines( Mat& img, const Point* const* pts, const int* npts, int ncontours, bool isClosed, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:function:: void polylines( InputOutputArray img, InputArrayOfArrays pts, bool isClosed, const Scalar& color, int thickness=1, int lineType=LINE_8, int shift=0 )
.. ocv:pyfunction:: cv2.polylines(img, pts, isClosed, color[, thickness[, lineType[, shift]]]) -> img
.. ocv:cfunction:: void cvPolyLine( CvArr* img, CvPoint** pts, const int* npts, int contours, int is_closed, CvScalar color, int thickness=1, int line_type=8, int shift=0 )
:param img: Image.
:param pts: Array of polygonal curves.
:param npts: Array of polygon vertex counters.
:param ncontours: Number of curves.
:param isClosed: Flag indicating whether the drawn polylines are closed or not. If they are closed, the function draws a line from the last vertex of each curve to its first vertex.
:param color: Polyline color.
:param thickness: Thickness of the polyline edges.
:param lineType: Type of the line segments. See the :ocv:func:`line` description.
:param shift: Number of fractional bits in the vertex coordinates.
The function ``polylines`` draws one or more polygonal curves.
drawContours
----------------
Draws contours outlines or filled contours.
.. ocv:function:: void drawContours( InputOutputArray image, InputArrayOfArrays contours, int contourIdx, const Scalar& color, int thickness=1, int lineType=LINE_8, InputArray hierarchy=noArray(), int maxLevel=INT_MAX, Point offset=Point() )
.. ocv:pyfunction:: cv2.drawContours(image, contours, contourIdx, color[, thickness[, lineType[, hierarchy[, maxLevel[, offset]]]]]) -> image
.. ocv:cfunction:: void cvDrawContours( CvArr * img, CvSeq* contour, CvScalar external_color, CvScalar hole_color, int max_level, int thickness=1, int line_type=8, CvPoint offset=cvPoint(0,0) )
:param image: Destination image.
:param contours: All the input contours. Each contour is stored as a point vector.
:param contourIdx: Parameter indicating a contour to draw. If it is negative, all the contours are drawn.
:param color: Color of the contours.
:param thickness: Thickness of lines the contours are drawn with. If it is negative (for example, ``thickness=CV_FILLED`` ), the contour interiors are
drawn.
:param lineType: Line connectivity. See :ocv:func:`line` for details.
:param hierarchy: Optional information about hierarchy. It is only needed if you want to draw only some of the contours (see ``maxLevel`` ).
:param maxLevel: Maximal level for drawn contours. If it is 0, only
the specified contour is drawn. If it is 1, the function draws the contour(s) and all the nested contours. If it is 2, the function draws the contours, all the nested contours, all the nested-to-nested contours, and so on. This parameter is only taken into account when there is ``hierarchy`` available.
:param offset: Optional contour shift parameter. Shift all the drawn contours by the specified :math:`\texttt{offset}=(dx,dy)` .
:param contour: Pointer to the first contour.
:param external_color: Color of external contours.
:param hole_color: Color of internal contours (holes).
The function draws contour outlines in the image if
:math:`\texttt{thickness} \ge 0` or fills the area bounded by the contours if
:math:`\texttt{thickness}<0` . The example below shows how to retrieve connected components from the binary image and label them: ::
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
using namespace cv;
using namespace std;
int main( int argc, char** argv )
{
Mat src;
// the first command-line parameter must be a filename of the binary
// (black-n-white) image
if( argc != 2 || !(src=imread(argv[1], 0)).data)
return -1;
Mat dst = Mat::zeros(src.rows, src.cols, CV_8UC3);
src = src > 1;
namedWindow( "Source", 1 );
imshow( "Source", src );
vector<vector<Point> > contours;
vector<Vec4i> hierarchy;
findContours( src, contours, hierarchy,
RETR_CCOMP, CHAIN_APPROX_SIMPLE );
// iterate through all the top-level contours,
// draw each connected component with its own random color
int idx = 0;
for( ; idx >= 0; idx = hierarchy[idx][0] )
{
Scalar color( rand()&255, rand()&255, rand()&255 );
drawContours( dst, contours, idx, color, FILLED, 8, hierarchy );
}
namedWindow( "Components", 1 );
imshow( "Components", dst );
waitKey(0);
}
.. note::
* An example using the drawContour functionality can be found at opencv_source_code/samples/cpp/contours2.cpp
* An example using drawContours to clean up a background segmentation result at opencv_source_code/samples/cpp/segment_objects.cpp
* (Python) An example using the drawContour functionality can be found at opencv_source/samples/python2/contours.py
putText
-----------
Draws a text string.
.. ocv:function:: void putText( InputOutputArray img, const String& text, Point org, int fontFace, double fontScale, Scalar color, int thickness=1, int lineType=LINE_8, bool bottomLeftOrigin=false )
.. ocv:pyfunction:: cv2.putText(img, text, org, fontFace, fontScale, color[, thickness[, lineType[, bottomLeftOrigin]]]) -> None
.. ocv:cfunction:: void cvPutText( CvArr* img, const char* text, CvPoint org, const CvFont* font, CvScalar color )
:param img: Image.
:param text: Text string to be drawn.
:param org: Bottom-left corner of the text string in the image.
:param font: ``CvFont`` structure initialized using :ocv:cfunc:`InitFont`.
:param fontFace: Font type. One of ``FONT_HERSHEY_SIMPLEX``, ``FONT_HERSHEY_PLAIN``, ``FONT_HERSHEY_DUPLEX``, ``FONT_HERSHEY_COMPLEX``, ``FONT_HERSHEY_TRIPLEX``, ``FONT_HERSHEY_COMPLEX_SMALL``, ``FONT_HERSHEY_SCRIPT_SIMPLEX``, or ``FONT_HERSHEY_SCRIPT_COMPLEX``,
where each of the font ID's can be combined with ``FONT_ITALIC`` to get the slanted letters.
:param fontScale: Font scale factor that is multiplied by the font-specific base size.
:param color: Text color.
:param thickness: Thickness of the lines used to draw a text.
:param lineType: Line type. See the ``line`` for details.
:param bottomLeftOrigin: When true, the image data origin is at the bottom-left corner. Otherwise, it is at the top-left corner.
The function ``putText`` renders the specified text string in the image.
Symbols that cannot be rendered using the specified font are
replaced by question marks. See
:ocv:func:`getTextSize` for a text rendering code example.

341
modules/core/doc/optim.rst Normal file
View File

@@ -0,0 +1,341 @@
Optimization Algorithms
=======================
.. highlight:: cpp
The algorithms in this section minimize or maximize function value within specified constraints or without any constraints.
solveLP
--------------------
Solve given (non-integer) linear programming problem using the Simplex Algorithm (Simplex Method).
What we mean here by "linear programming problem" (or LP problem, for short) can be
formulated as:
.. math::
\mbox{Maximize } c\cdot x\\
\mbox{Subject to:}\\
Ax\leq b\\
x\geq 0
Where :math:`c` is fixed *1*-by-*n* row-vector, :math:`A` is fixed *m*-by-*n* matrix, :math:`b` is fixed *m*-by-*1* column vector and
:math:`x` is an arbitrary *n*-by-*1* column vector, which satisfies the constraints.
Simplex algorithm is one of many algorithms that are designed to handle this sort of problems efficiently. Although it is not optimal in theoretical
sense (there exist algorithms that can solve any problem written as above in polynomial type, while simplex method degenerates to exponential time
for some special cases), it is well-studied, easy to implement and is shown to work well for real-life purposes.
The particular implementation is taken almost verbatim from **Introduction to Algorithms, third edition**
by T. H. Cormen, C. E. Leiserson, R. L. Rivest and Clifford Stein. In particular, the Bland's rule
(`http://en.wikipedia.org/wiki/Bland%27s\_rule <http://en.wikipedia.org/wiki/Bland%27s_rule>`_) is used to prevent cycling.
.. ocv:function:: int solveLP(const Mat& Func, const Mat& Constr, Mat& z)
:param Func: This row-vector corresponds to :math:`c` in the LP problem formulation (see above). It should contain 32- or 64-bit floating point numbers. As a convenience, column-vector may be also submitted, in the latter case it is understood to correspond to :math:`c^T`.
:param Constr: *m*-by-*n\+1* matrix, whose rightmost column corresponds to :math:`b` in formulation above and the remaining to :math:`A`. It should containt 32- or 64-bit floating point numbers.
:param z: The solution will be returned here as a column-vector - it corresponds to :math:`c` in the formulation above. It will contain 64-bit floating point numbers.
:return: One of the return codes:
::
//!the return codes for solveLP() function
enum
{
SOLVELP_UNBOUNDED = -2, //problem is unbounded (target function can achieve arbitrary high values)
SOLVELP_UNFEASIBLE = -1, //problem is unfeasible (there are no points that satisfy all the constraints imposed)
SOLVELP_SINGLE = 0, //there is only one maximum for target function
SOLVELP_MULTI = 1 //there are multiple maxima for target function - the arbitrary one is returned
};
DownhillSolver
---------------------------------
.. ocv:class:: DownhillSolver
This class is used to perform the non-linear non-constrained *minimization* of a function, defined on an *n*-dimensional Euclidean space,
using the **Nelder-Mead method**, also known as **downhill simplex method**. The basic idea about the method can be obtained from
(`http://en.wikipedia.org/wiki/Nelder-Mead\_method <http://en.wikipedia.org/wiki/Nelder-Mead_method>`_). It should be noted, that
this method, although deterministic, is rather a heuristic and therefore may converge to a local minima, not necessary a global one.
It is iterative optimization technique, which at each step uses an information about the values of a function evaluated only at
*n+1* points, arranged as a *simplex* in *n*-dimensional space (hence the second name of the method). At each step new point is
chosen to evaluate function at, obtained value is compared with previous ones and based on this information simplex changes it's shape
, slowly moving to the local minimum. Thus this method is using *only* function values to make decision, on contrary to, say, Nonlinear
Conjugate Gradient method (which is also implemented in ``optim``).
Algorithm stops when the number of function evaluations done exceeds ``termcrit.maxCount``, when the function values at the
vertices of simplex are within ``termcrit.epsilon`` range or simplex becomes so small that it
can enclosed in a box with ``termcrit.epsilon`` sides, whatever comes first, for some defined by user
positive integer ``termcrit.maxCount`` and positive non-integer ``termcrit.epsilon``.
::
class CV_EXPORTS Solver : public Algorithm
{
public:
class CV_EXPORTS Function
{
public:
virtual ~Function() {}
virtual double calc(const double* x) const = 0;
virtual void getGradient(const double* /*x*/,double* /*grad*/) {}
};
virtual Ptr<Function> getFunction() const = 0;
virtual void setFunction(const Ptr<Function>& f) = 0;
virtual TermCriteria getTermCriteria() const = 0;
virtual void setTermCriteria(const TermCriteria& termcrit) = 0;
// x contain the initial point before the call and the minima position (if algorithm converged) after. x is assumed to be (something that
// after getMat() will return) row-vector or column-vector. *It's size and should
// be consisted with previous dimensionality data given, if any (otherwise, it determines dimensionality)*
virtual double minimize(InputOutputArray x) = 0;
};
class CV_EXPORTS DownhillSolver : public Solver
{
public:
//! returns row-vector, even if the column-vector was given
virtual void getInitStep(OutputArray step) const=0;
//!This should be called at least once before the first call to minimize() and step is assumed to be (something that
//! after getMat() will return) row-vector or column-vector. *It's dimensionality determines the dimensionality of a problem.*
virtual void setInitStep(InputArray step)=0;
};
It should be noted, that ``DownhillSolver`` is a derivative of the abstract interface ``Solver``, which in
turn is derived from the ``Algorithm`` interface and is used to encapsulate the functionality, common to all non-linear optimization
algorithms in the ``optim`` module.
DownhillSolver::getFunction
--------------------------------------------
Getter for the optimized function. The optimized function is represented by ``Solver::Function`` interface, which requires
derivatives to implement the sole method ``calc(double*)`` to evaluate the function.
.. ocv:function:: Ptr<Solver::Function> DownhillSolver::getFunction()
:return: Smart-pointer to an object that implements ``Solver::Function`` interface - it represents the function that is being optimized. It can be empty, if no function was given so far.
DownhillSolver::setFunction
-----------------------------------------------
Setter for the optimized function. *It should be called at least once before the call to* ``DownhillSolver::minimize()``, as
default value is not usable.
.. ocv:function:: void DownhillSolver::setFunction(const Ptr<Solver::Function>& f)
:param f: The new function to optimize.
DownhillSolver::getTermCriteria
----------------------------------------------------
Getter for the previously set terminal criteria for this algorithm.
.. ocv:function:: TermCriteria DownhillSolver::getTermCriteria()
:return: Deep copy of the terminal criteria used at the moment.
DownhillSolver::setTermCriteria
------------------------------------------
Set terminal criteria for downhill simplex method. Two things should be noted. First, this method *is not necessary* to be called
before the first call to ``DownhillSolver::minimize()``, as the default value is sensible. Second, the method will raise an error
if ``termcrit.type!=(TermCriteria::MAX_ITER+TermCriteria::EPS)``, ``termcrit.epsilon<=0`` or ``termcrit.maxCount<=0``. That is,
both ``epsilon`` and ``maxCount`` should be set to positive values (non-integer and integer respectively) and they represent
tolerance and maximal number of function evaluations that is allowed.
Algorithm stops when the number of function evaluations done exceeds ``termcrit.maxCount``, when the function values at the
vertices of simplex are within ``termcrit.epsilon`` range or simplex becomes so small that it
can enclosed in a box with ``termcrit.epsilon`` sides, whatever comes first.
.. ocv:function:: void DownhillSolver::setTermCriteria(const TermCriteria& termcrit)
:param termcrit: Terminal criteria to be used, represented as ``TermCriteria`` structure (defined elsewhere in openCV). Mind you, that it should meet ``(termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0)``, otherwise the error will be raised.
DownhillSolver::getInitStep
-----------------------------------
Returns the initial step that will be used in downhill simplex algorithm. See the description
of corresponding setter (follows next) for the meaning of this parameter.
.. ocv:function:: void getInitStep(OutputArray step)
:param step: Initial step that will be used in algorithm. Note, that although corresponding setter accepts column-vectors as well as row-vectors, this method will return a row-vector.
DownhillSolver::setInitStep
----------------------------------
Sets the initial step that will be used in downhill simplex algorithm. Step, together with initial point (givin in ``DownhillSolver::minimize``)
are two *n*-dimensional vectors that are used to determine the shape of initial simplex. Roughly said, initial point determines the position
of a simplex (it will become simplex's centroid), while step determines the spread (size in each dimension) of a simplex. To be more precise,
if :math:`s,x_0\in\mathbb{R}^n` are the initial step and initial point respectively, the vertices of a simplex will be: :math:`v_0:=x_0-\frac{1}{2}
s` and :math:`v_i:=x_0+s_i` for :math:`i=1,2,\dots,n` where :math:`s_i` denotes projections of the initial step of *n*-th coordinate (the result
of projection is treated to be vector given by :math:`s_i:=e_i\cdot\left<e_i\cdot s\right>`, where :math:`e_i` form canonical basis)
.. ocv:function:: void setInitStep(InputArray step)
:param step: Initial step that will be used in algorithm. Roughly said, it determines the spread (size in each dimension) of an initial simplex.
DownhillSolver::minimize
-----------------------------------
The main method of the ``DownhillSolver``. It actually runs the algorithm and performs the minimization. The sole input parameter determines the
centroid of the starting simplex (roughly, it tells where to start), all the others (terminal criteria, initial step, function to be minimized)
are supposed to be set via the setters before the call to this method or the default values (not always sensible) will be used.
.. ocv:function:: double DownhillSolver::minimize(InputOutputArray x)
:param x: The initial point, that will become a centroid of an initial simplex. After the algorithm will terminate, it will be setted to the point where the algorithm stops, the point of possible minimum.
:return: The value of a function at the point found.
createDownhillSolver
------------------------------------
This function returns the reference to the ready-to-use ``DownhillSolver`` object. All the parameters are optional, so this procedure can be called
even without parameters at all. In this case, the default values will be used. As default value for terminal criteria are the only sensible ones,
``DownhillSolver::setFunction()`` and ``DownhillSolver::setInitStep()`` should be called upon the obtained object, if the respective parameters
were not given to ``createDownhillSolver()``. Otherwise, the two ways (give parameters to ``createDownhillSolver()`` or miss them out and call the
``DownhillSolver::setFunction()`` and ``DownhillSolver::setInitStep()``) are absolutely equivalent (and will drop the same errors in the same way,
should invalid input be detected).
.. ocv:function:: Ptr<DownhillSolver> createDownhillSolver(const Ptr<Solver::Function>& f,InputArray initStep, TermCriteria termcrit)
:param f: Pointer to the function that will be minimized, similarly to the one you submit via ``DownhillSolver::setFunction``.
:param step: Initial step, that will be used to construct the initial simplex, similarly to the one you submit via ``DownhillSolver::setInitStep``.
:param termcrit: Terminal criteria to the algorithm, similarly to the one you submit via ``DownhillSolver::setTermCriteria``.
ConjGradSolver
---------------------------------
.. ocv:class:: ConjGradSolver
This class is used to perform the non-linear non-constrained *minimization* of a function with *known gradient*
, defined on an *n*-dimensional Euclidean space,
using the **Nonlinear Conjugate Gradient method**. The implementation was done based on the beautifully clear explanatory article `An Introduction to the Conjugate Gradient Method Without the Agonizing Pain <http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf>`_
by Jonathan Richard Shewchuk. The method can be seen as an adaptation of a standard Conjugate Gradient method (see, for example
`http://en.wikipedia.org/wiki/Conjugate_gradient_method <http://en.wikipedia.org/wiki/Conjugate_gradient_method>`_) for numerically solving the
systems of linear equations.
It should be noted, that
this method, although deterministic, is rather a heuristic method and therefore may converge to a local minima, not necessary a global one. What
is even more disastrous, most of its behaviour is ruled by gradient, therefore it essentially cannot distinguish between local minima and maxima.
Therefore, if it starts sufficiently near to the local maximum, it may converge to it. Another obvious restriction is that it should be possible
to compute the gradient of a function at any point, thus it is preferable to have analytic expression for gradient and computational burden
should be born by the user.
The latter responsibility is accompilished via the ``getGradient(const double* x,double* grad)`` method of a
``Solver::Function`` interface (which represents function that is being optimized). This method takes point a point in *n*-dimensional space
(first argument represents the array of coordinates of that point) and comput its gradient (it should be stored in the second argument as an array).
::
class CV_EXPORTS Solver : public Algorithm
{
public:
class CV_EXPORTS Function
{
public:
virtual ~Function() {}
virtual double calc(const double* x) const = 0;
virtual void getGradient(const double* /*x*/,double* /*grad*/) {}
};
virtual Ptr<Function> getFunction() const = 0;
virtual void setFunction(const Ptr<Function>& f) = 0;
virtual TermCriteria getTermCriteria() const = 0;
virtual void setTermCriteria(const TermCriteria& termcrit) = 0;
// x contain the initial point before the call and the minima position (if algorithm converged) after. x is assumed to be (something that
// after getMat() will return) row-vector or column-vector. *It's size and should
// be consisted with previous dimensionality data given, if any (otherwise, it determines dimensionality)*
virtual double minimize(InputOutputArray x) = 0;
};
class CV_EXPORTS ConjGradSolver : public Solver{
};
Note, that class ``ConjGradSolver`` thus does not add any new methods to the basic ``Solver`` interface.
ConjGradSolver::getFunction
--------------------------------------------
Getter for the optimized function. The optimized function is represented by ``Solver::Function`` interface, which requires
derivatives to implement the method ``calc(double*)`` to evaluate the function. It should be emphasized once more, that since Nonlinear
Conjugate Gradient method requires gradient to be computable in addition to the function values,
``getGradient(const double* x,double* grad)`` method of a ``Solver::Function`` interface should be also implemented meaningfully.
.. ocv:function:: Ptr<Solver::Function> ConjGradSolver::getFunction()
:return: Smart-pointer to an object that implements ``Solver::Function`` interface - it represents the function that is being optimized. It can be empty, if no function was given so far.
ConjGradSolver::setFunction
-----------------------------------------------
Setter for the optimized function. *It should be called at least once before the call to* ``ConjGradSolver::minimize()``, as
default value is not usable.
.. ocv:function:: void ConjGradSolver::setFunction(const Ptr<Solver::Function>& f)
:param f: The new function to optimize.
ConjGradSolver::getTermCriteria
----------------------------------------------------
Getter for the previously set terminal criteria for this algorithm.
.. ocv:function:: TermCriteria ConjGradSolver::getTermCriteria()
:return: Deep copy of the terminal criteria used at the moment.
ConjGradSolver::setTermCriteria
------------------------------------------
Set terminal criteria for downhill simplex method. Two things should be noted. First, this method *is not necessary* to be called
before the first call to ``ConjGradSolver::minimize()``, as the default value is sensible. Second, the method will raise an error
if ``termcrit.type!=(TermCriteria::MAX_ITER+TermCriteria::EPS)`` and ``termcrit.type!=TermCriteria::MAX_ITER``. This means that termination criteria
has to restrict maximum number of iterations to be done and may optionally allow algorithm to stop earlier if certain tolerance
is achieved (what we mean by "tolerance is achieved" will be clarified below). If ``termcrit`` restricts both tolerance and maximum iteration
number, both ``termcrit.epsilon`` and ``termcrit.maxCount`` should be positive. In case, if ``termcrit.type==TermCriteria::MAX_ITER``,
only member ``termcrit.maxCount`` is required to be positive and in this case algorithm will just work for required number of iterations.
In current implementation, "tolerance is achieved" means that we have arrived at the point where the :math:`L_2`-norm of the gradient is less
than the tolerance value.
.. ocv:function:: void ConjGradSolver::setTermCriteria(const TermCriteria& termcrit)
:param termcrit: Terminal criteria to be used, represented as ``TermCriteria`` structure (defined elsewhere in openCV). Mind you, that it should meet ``termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0`` or ``termcrit.type==TermCriteria::MAX_ITER) && termcrit.maxCount>0``, otherwise the error will be raised.
ConjGradSolver::minimize
-----------------------------------
The main method of the ``ConjGradSolver``. It actually runs the algorithm and performs the minimization. The sole input parameter determines the
centroid of the starting simplex (roughly, it tells where to start), all the others (terminal criteria and function to be minimized)
are supposed to be set via the setters before the call to this method or the default values (not always sensible) will be used. Sometimes it may
throw an error, if these default values cannot be used (say, you forgot to set the function to minimize and default value, that is, empty function,
cannot be used).
.. ocv:function:: double ConjGradSolver::minimize(InputOutputArray x)
:param x: The initial point. It is hard to overemphasize how important the choise of initial point is when you are using the heuristic algorithm like this one. Badly chosen initial point can make algorithm converge to (local) maximum instead of minimum, do not converge at all, converge to local minimum instead of global one.
:return: The value of a function at the point found.
createConjGradSolver
------------------------------------
This function returns the reference to the ready-to-use ``ConjGradSolver`` object. All the parameters are optional, so this procedure can be called
even without parameters at all. In this case, the default values will be used. As default value for terminal criteria are the only sensible ones,
``ConjGradSolver::setFunction()`` should be called upon the obtained object, if the function
was not given to ``createConjGradSolver()``. Otherwise, the two ways (submit it to ``createConjGradSolver()`` or miss it out and call the
``ConjGradSolver::setFunction()``) are absolutely equivalent (and will drop the same errors in the same way,
should invalid input be detected).
.. ocv:function:: Ptr<ConjGradSolver> createConjGradSolver(const Ptr<Solver::Function>& f, TermCriteria termcrit)
:param f: Pointer to the function that will be minimized, similarly to the one you submit via ``ConjGradSolver::setFunction``.
:param termcrit: Terminal criteria to the algorithm, similarly to the one you submit via ``ConjGradSolver::setTermCriteria``.

BIN
modules/core/doc/pics/ellipse.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 2.4 KiB

92
modules/core/include/opencv2/core.hpp
View File

@@ -506,96 +506,6 @@ CV_EXPORTS_W void randn(InputOutputArray dst, InputArray mean, InputArray stddev
//! shuffles the input array elements
CV_EXPORTS_W void randShuffle(InputOutputArray dst, double iterFactor = 1., RNG* rng = 0);
//! draws the line segment (pt1, pt2) in the image
CV_EXPORTS_W void line(InputOutputArray img, Point pt1, Point pt2, const Scalar& color,
int thickness = 1, int lineType = LINE_8, int shift = 0);
//! draws an arrow from pt1 to pt2 in the image
CV_EXPORTS_W void arrowedLine(InputOutputArray img, Point pt1, Point pt2, const Scalar& color,
int thickness=1, int line_type=8, int shift=0, double tipLength=0.1);
//! draws the rectangle outline or a solid rectangle with the opposite corners pt1 and pt2 in the image
CV_EXPORTS_W void rectangle(InputOutputArray img, Point pt1, Point pt2,
const Scalar& color, int thickness = 1,
int lineType = LINE_8, int shift = 0);
//! draws the rectangle outline or a solid rectangle covering rec in the image
CV_EXPORTS void rectangle(CV_IN_OUT Mat& img, Rect rec,
const Scalar& color, int thickness = 1,
int lineType = LINE_8, int shift = 0);
//! draws the circle outline or a solid circle in the image
CV_EXPORTS_W void circle(InputOutputArray img, Point center, int radius,
const Scalar& color, int thickness = 1,
int lineType = LINE_8, int shift = 0);
//! draws an elliptic arc, ellipse sector or a rotated ellipse in the image
CV_EXPORTS_W void ellipse(InputOutputArray img, Point center, Size axes,
double angle, double startAngle, double endAngle,
const Scalar& color, int thickness = 1,
int lineType = LINE_8, int shift = 0);
//! draws a rotated ellipse in the image
CV_EXPORTS_W void ellipse(InputOutputArray img, const RotatedRect& box, const Scalar& color,
int thickness = 1, int lineType = LINE_8);
//! draws a filled convex polygon in the image
CV_EXPORTS void fillConvexPoly(Mat& img, const Point* pts, int npts,
const Scalar& color, int lineType = LINE_8,
int shift = 0);
CV_EXPORTS_W void fillConvexPoly(InputOutputArray img, InputArray points,
const Scalar& color, int lineType = LINE_8,
int shift = 0);
//! fills an area bounded by one or more polygons
CV_EXPORTS void fillPoly(Mat& img, const Point** pts,
const int* npts, int ncontours,
const Scalar& color, int lineType = LINE_8, int shift = 0,
Point offset = Point() );
CV_EXPORTS_W void fillPoly(InputOutputArray img, InputArrayOfArrays pts,
const Scalar& color, int lineType = LINE_8, int shift = 0,
Point offset = Point() );
//! draws one or more polygonal curves
CV_EXPORTS void polylines(Mat& img, const Point* const* pts, const int* npts,
int ncontours, bool isClosed, const Scalar& color,
int thickness = 1, int lineType = LINE_8, int shift = 0 );
CV_EXPORTS_W void polylines(InputOutputArray img, InputArrayOfArrays pts,
bool isClosed, const Scalar& color,
int thickness = 1, int lineType = LINE_8, int shift = 0 );
//! draws contours in the image
CV_EXPORTS_W void drawContours( InputOutputArray image, InputArrayOfArrays contours,
int contourIdx, const Scalar& color,
int thickness = 1, int lineType = LINE_8,
InputArray hierarchy = noArray(),
int maxLevel = INT_MAX, Point offset = Point() );
//! clips the line segment by the rectangle Rect(0, 0, imgSize.width, imgSize.height)
CV_EXPORTS bool clipLine(Size imgSize, CV_IN_OUT Point& pt1, CV_IN_OUT Point& pt2);
//! clips the line segment by the rectangle imgRect
CV_EXPORTS_W bool clipLine(Rect imgRect, CV_OUT CV_IN_OUT Point& pt1, CV_OUT CV_IN_OUT Point& pt2);
//! converts elliptic arc to a polygonal curve
CV_EXPORTS_W void ellipse2Poly( Point center, Size axes, int angle,
int arcStart, int arcEnd, int delta,
CV_OUT std::vector<Point>& pts );
//! renders text string in the image
CV_EXPORTS_W void putText( InputOutputArray img, const String& text, Point org,
int fontFace, double fontScale, Scalar color,
int thickness = 1, int lineType = LINE_8,
bool bottomLeftOrigin = false );
//! returns bounding box of the text string
CV_EXPORTS_W Size getTextSize(const String& text, int fontFace,
double fontScale, int thickness,
CV_OUT int* baseLine);
/*!
Principal Component Analysis
@@ -1319,5 +1229,7 @@ template<> struct ParamType<uchar>
#include "opencv2/core/operations.hpp"
#include "opencv2/core/cvstd.inl.hpp"
#include "opencv2/core/utility.hpp"
#include "opencv2/core/optim.hpp"
#endif /*__OPENCV_CORE_HPP__*/

186
modules/core/include/opencv2/core/core_c.h
View File

@@ -1262,192 +1262,6 @@ CVAPI(int) cvNextGraphItem( CvGraphScanner* scanner );
/* Creates a copy of graph */
CVAPI(CvGraph*) cvCloneGraph( const CvGraph* graph, CvMemStorage* storage );
/****************************************************************************************\
* Drawing *
\****************************************************************************************/
/****************************************************************************************\
* Drawing functions work with images/matrices of arbitrary type. *
* For color images the channel order is BGR[A] *
* Antialiasing is supported only for 8-bit image now. *
* All the functions include parameter color that means rgb value (that may be *
* constructed with CV_RGB macro) for color images and brightness *
* for grayscale images. *
* If a drawn figure is partially or completely outside of the image, it is clipped.*
\****************************************************************************************/
#define CV_RGB( r, g, b ) cvScalar( (b), (g), (r), 0 )
#define CV_FILLED -1
#define CV_AA 16
/* Draws 4-connected, 8-connected or antialiased line segment connecting two points */
CVAPI(void) cvLine( CvArr* img, CvPoint pt1, CvPoint pt2,
CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0) );
/* Draws a rectangle given two opposite corners of the rectangle (pt1 & pt2),
if thickness<0 (e.g. thickness == CV_FILLED), the filled box is drawn */
CVAPI(void) cvRectangle( CvArr* img, CvPoint pt1, CvPoint pt2,
CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8),
int shift CV_DEFAULT(0));
/* Draws a rectangle specified by a CvRect structure */
CVAPI(void) cvRectangleR( CvArr* img, CvRect r,
CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8),
int shift CV_DEFAULT(0));
/* Draws a circle with specified center and radius.
Thickness works in the same way as with cvRectangle */
CVAPI(void) cvCircle( CvArr* img, CvPoint center, int radius,
CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0));
/* Draws ellipse outline, filled ellipse, elliptic arc or filled elliptic sector,
depending on <thickness>, <start_angle> and <end_angle> parameters. The resultant figure
is rotated by <angle>. All the angles are in degrees */
CVAPI(void) cvEllipse( CvArr* img, CvPoint center, CvSize axes,
double angle, double start_angle, double end_angle,
CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0));
CV_INLINE void cvEllipseBox( CvArr* img, CvBox2D box, CvScalar color,
int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0) )
{
CvSize axes;
axes.width = cvRound(box.size.width*0.5);
axes.height = cvRound(box.size.height*0.5);
cvEllipse( img, cvPointFrom32f( box.center ), axes, box.angle,
0, 360, color, thickness, line_type, shift );
}
/* Fills convex or monotonous polygon. */
CVAPI(void) cvFillConvexPoly( CvArr* img, const CvPoint* pts, int npts, CvScalar color,
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0));
/* Fills an area bounded by one or more arbitrary polygons */
CVAPI(void) cvFillPoly( CvArr* img, CvPoint** pts, const int* npts,
int contours, CvScalar color,
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0) );
/* Draws one or more polygonal curves */
CVAPI(void) cvPolyLine( CvArr* img, CvPoint** pts, const int* npts, int contours,
int is_closed, CvScalar color, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8), int shift CV_DEFAULT(0) );
#define cvDrawRect cvRectangle
#define cvDrawLine cvLine
#define cvDrawCircle cvCircle
#define cvDrawEllipse cvEllipse
#define cvDrawPolyLine cvPolyLine
/* Clips the line segment connecting *pt1 and *pt2
by the rectangular window
(0<=x<img_size.width, 0<=y<img_size.height). */
CVAPI(int) cvClipLine( CvSize img_size, CvPoint* pt1, CvPoint* pt2 );
/* Initializes line iterator. Initially, line_iterator->ptr will point
to pt1 (or pt2, see left_to_right description) location in the image.
Returns the number of pixels on the line between the ending points. */
CVAPI(int) cvInitLineIterator( const CvArr* image, CvPoint pt1, CvPoint pt2,
CvLineIterator* line_iterator,
int connectivity CV_DEFAULT(8),
int left_to_right CV_DEFAULT(0));
/* Moves iterator to the next line point */
#define CV_NEXT_LINE_POINT( line_iterator ) \
{ \
int _line_iterator_mask = (line_iterator).err < 0 ? -1 : 0; \
(line_iterator).err += (line_iterator).minus_delta + \
((line_iterator).plus_delta & _line_iterator_mask); \
(line_iterator).ptr += (line_iterator).minus_step + \
((line_iterator).plus_step & _line_iterator_mask); \
}
/* basic font types */
#define CV_FONT_HERSHEY_SIMPLEX 0
#define CV_FONT_HERSHEY_PLAIN 1
#define CV_FONT_HERSHEY_DUPLEX 2
#define CV_FONT_HERSHEY_COMPLEX 3
#define CV_FONT_HERSHEY_TRIPLEX 4
#define CV_FONT_HERSHEY_COMPLEX_SMALL 5
#define CV_FONT_HERSHEY_SCRIPT_SIMPLEX 6
#define CV_FONT_HERSHEY_SCRIPT_COMPLEX 7
/* font flags */
#define CV_FONT_ITALIC 16
#define CV_FONT_VECTOR0 CV_FONT_HERSHEY_SIMPLEX
/* Font structure */
typedef struct CvFont
{
const char* nameFont; //Qt:nameFont
CvScalar color; //Qt:ColorFont -> cvScalar(blue_component, green_component, red\_component[, alpha_component])
int font_face; //Qt: bool italic /* =CV_FONT_* */
const int* ascii; /* font data and metrics */
const int* greek;
const int* cyrillic;
float hscale, vscale;
float shear; /* slope coefficient: 0 - normal, >0 - italic */
int thickness; //Qt: weight /* letters thickness */
float dx; /* horizontal interval between letters */
int line_type; //Qt: PointSize
}
CvFont;
/* Initializes font structure used further in cvPutText */
CVAPI(void) cvInitFont( CvFont* font, int font_face,
double hscale, double vscale,
double shear CV_DEFAULT(0),
int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8));
CV_INLINE CvFont cvFont( double scale, int thickness CV_DEFAULT(1) )
{
CvFont font;
cvInitFont( &font, CV_FONT_HERSHEY_PLAIN, scale, scale, 0, thickness, CV_AA );
return font;
}
/* Renders text stroke with specified font and color at specified location.
CvFont should be initialized with cvInitFont */
CVAPI(void) cvPutText( CvArr* img, const char* text, CvPoint org,
const CvFont* font, CvScalar color );
/* Calculates bounding box of text stroke (useful for alignment) */
CVAPI(void) cvGetTextSize( const char* text_string, const CvFont* font,
CvSize* text_size, int* baseline );
/* Unpacks color value, if arrtype is CV_8UC?, <color> is treated as
packed color value, otherwise the first channels (depending on arrtype)
of destination scalar are set to the same value = <color> */
CVAPI(CvScalar) cvColorToScalar( double packed_color, int arrtype );
/* Returns the polygon points which make up the given ellipse. The ellipse is define by
the box of size 'axes' rotated 'angle' around the 'center'. A partial sweep
of the ellipse arc can be done by spcifying arc_start and arc_end to be something
other than 0 and 360, respectively. The input array 'pts' must be large enough to
hold the result. The total number of points stored into 'pts' is returned by this
function. */
CVAPI(int) cvEllipse2Poly( CvPoint center, CvSize axes,
int angle, int arc_start, int arc_end, CvPoint * pts, int delta );
/* Draws contour outlines or filled interiors on the image */
CVAPI(void) cvDrawContours( CvArr *img, CvSeq* contour,
CvScalar external_color, CvScalar hole_color,
int max_level, int thickness CV_DEFAULT(1),
int line_type CV_DEFAULT(8),
CvPoint offset CV_DEFAULT(cvPoint(0,0)));
/* Does look-up transformation. Elements of the source array
(that should be 8uC1 or 8sC1) are used as indexes in lutarr 256-element table */

111
modules/core/include/opencv2/core/optim.hpp Normal file
View File

@@ -0,0 +1,111 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#ifndef __OPENCV_OPTIM_HPP__
#define __OPENCV_OPTIM_HPP__
#include "opencv2/core.hpp"
namespace cv
{
class CV_EXPORTS MinProblemSolver : public Algorithm
{
public:
class CV_EXPORTS Function
{
public:
virtual ~Function() {}
virtual double calc(const double* x) const = 0;
virtual void getGradient(const double* /*x*/,double* /*grad*/) {}
};
virtual Ptr<Function> getFunction() const = 0;
virtual void setFunction(const Ptr<Function>& f) = 0;
virtual TermCriteria getTermCriteria() const = 0;
virtual void setTermCriteria(const TermCriteria& termcrit) = 0;
// x contain the initial point before the call and the minima position (if algorithm converged) after. x is assumed to be (something that
// after getMat() will return) row-vector or column-vector. *It's size and should
// be consisted with previous dimensionality data given, if any (otherwise, it determines dimensionality)*
virtual double minimize(InputOutputArray x) = 0;
};
//! downhill simplex class
class CV_EXPORTS DownhillSolver : public MinProblemSolver
{
public:
//! returns row-vector, even if the column-vector was given
virtual void getInitStep(OutputArray step) const=0;
//!This should be called at least once before the first call to minimize() and step is assumed to be (something that
//! after getMat() will return) row-vector or column-vector. *It's dimensionality determines the dimensionality of a problem.*
virtual void setInitStep(InputArray step)=0;
// both minRange & minError are specified by termcrit.epsilon;
// In addition, user may specify the number of iterations that the algorithm does.
static Ptr<DownhillSolver> create(const Ptr<MinProblemSolver::Function>& f=Ptr<MinProblemSolver::Function>(),
InputArray initStep=Mat_<double>(1,1,0.0),
TermCriteria termcrit=TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS,5000,0.000001));
};
//! conjugate gradient method
class CV_EXPORTS ConjGradSolver : public MinProblemSolver
{
public:
static Ptr<ConjGradSolver> create(const Ptr<MinProblemSolver::Function>& f=Ptr<ConjGradSolver::Function>(),
TermCriteria termcrit=TermCriteria(TermCriteria::MAX_ITER+TermCriteria::EPS,5000,0.000001));
};
//!the return codes for solveLP() function
enum
{
SOLVELP_UNBOUNDED = -2, //problem is unbounded (target function can achieve arbitrary high values)
SOLVELP_UNFEASIBLE = -1, //problem is unfeasible (there are no points that satisfy all the constraints imposed)
SOLVELP_SINGLE = 0, //there is only one maximum for target function
SOLVELP_MULTI = 1 //there are multiple maxima for target function - the arbitrary one is returned
};
CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z);
}// cv
#endif

186
modules/core/src/conjugate_gradient.cpp Normal file
View File

@@ -0,0 +1,186 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#define dprintf(x)
#define print_matrix(x)
namespace cv
{
#define SEC_METHOD_ITERATIONS 4
#define INITIAL_SEC_METHOD_SIGMA 0.1
class ConjGradSolverImpl : public ConjGradSolver
{
public:
Ptr<Function> getFunction() const;
void setFunction(const Ptr<Function>& f);
TermCriteria getTermCriteria() const;
ConjGradSolverImpl();
void setTermCriteria(const TermCriteria& termcrit);
double minimize(InputOutputArray x);
protected:
Ptr<MinProblemSolver::Function> _Function;
TermCriteria _termcrit;
Mat_<double> d,r,buf_x,r_old;
Mat_<double> minimizeOnTheLine_buf1,minimizeOnTheLine_buf2;
private:
static void minimizeOnTheLine(Ptr<MinProblemSolver::Function> _f,Mat_<double>& x,const Mat_<double>& d,Mat_<double>& buf1,Mat_<double>& buf2);
};
void ConjGradSolverImpl::minimizeOnTheLine(Ptr<MinProblemSolver::Function> _f,Mat_<double>& x,const Mat_<double>& d,Mat_<double>& buf1,
Mat_<double>& buf2){
double sigma=INITIAL_SEC_METHOD_SIGMA;
buf1=0.0;
buf2=0.0;
dprintf(("before minimizeOnTheLine\n"));
dprintf(("x:\n"));
print_matrix(x);
dprintf(("d:\n"));
print_matrix(d);
for(int i=0;i<SEC_METHOD_ITERATIONS;i++){
_f->getGradient((double*)x.data,(double*)buf1.data);
dprintf(("buf1:\n"));
print_matrix(buf1);
x=x+sigma*d;
_f->getGradient((double*)x.data,(double*)buf2.data);
dprintf(("buf2:\n"));
print_matrix(buf2);
double d1=buf1.dot(d), d2=buf2.dot(d);
if((d1-d2)==0){
break;
}
double alpha=-sigma*d1/(d2-d1);
dprintf(("(buf2.dot(d)-buf1.dot(d))=%f\nalpha=%f\n",(buf2.dot(d)-buf1.dot(d)),alpha));
x=x+(alpha-sigma)*d;
sigma=-alpha;
}
dprintf(("after minimizeOnTheLine\n"));
print_matrix(x);
}
double ConjGradSolverImpl::minimize(InputOutputArray x){
CV_Assert(_Function.empty()==false);
dprintf(("termcrit:\n\ttype: %d\n\tmaxCount: %d\n\tEPS: %g\n",_termcrit.type,_termcrit.maxCount,_termcrit.epsilon));
Mat x_mat=x.getMat();
CV_Assert(MIN(x_mat.rows,x_mat.cols)==1);
int ndim=MAX(x_mat.rows,x_mat.cols);
CV_Assert(x_mat.type()==CV_64FC1);
if(d.cols!=ndim){
d.create(1,ndim);
r.create(1,ndim);
r_old.create(1,ndim);
minimizeOnTheLine_buf1.create(1,ndim);
minimizeOnTheLine_buf2.create(1,ndim);
}
Mat_<double> proxy_x;
if(x_mat.rows>1){
buf_x.create(1,ndim);
Mat_<double> proxy(ndim,1,buf_x.ptr<double>());
x_mat.copyTo(proxy);
proxy_x=buf_x;
}else{
proxy_x=x_mat;
}
_Function->getGradient(proxy_x.ptr<double>(),d.ptr<double>());
d*=-1.0;
d.copyTo(r);
//here everything goes. check that everything is setted properly
dprintf(("proxy_x\n"));print_matrix(proxy_x);
dprintf(("d first time\n"));print_matrix(d);
dprintf(("r\n"));print_matrix(r);
double beta=0;
for(int count=0;count<_termcrit.maxCount;count++){
minimizeOnTheLine(_Function,proxy_x,d,minimizeOnTheLine_buf1,minimizeOnTheLine_buf2);
r.copyTo(r_old);
_Function->getGradient(proxy_x.ptr<double>(),r.ptr<double>());
r*=-1.0;
double r_norm_sq=norm(r);
if(_termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && r_norm_sq<_termcrit.epsilon){
break;
}
r_norm_sq=r_norm_sq*r_norm_sq;
beta=MAX(0.0,(r_norm_sq-r.dot(r_old))/r_norm_sq);
d=r+beta*d;
}
if(x_mat.rows>1){
Mat(ndim, 1, CV_64F, proxy_x.ptr<double>()).copyTo(x);
}
return _Function->calc(proxy_x.ptr<double>());
}
ConjGradSolverImpl::ConjGradSolverImpl(){
_Function=Ptr<Function>();
}
Ptr<MinProblemSolver::Function> ConjGradSolverImpl::getFunction()const{
return _Function;
}
void ConjGradSolverImpl::setFunction(const Ptr<Function>& f){
_Function=f;
}
TermCriteria ConjGradSolverImpl::getTermCriteria()const{
return _termcrit;
}
void ConjGradSolverImpl::setTermCriteria(const TermCriteria& termcrit){
CV_Assert((termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0) ||
((termcrit.type==TermCriteria::MAX_ITER) && termcrit.maxCount>0));
_termcrit=termcrit;
}
// both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does.
Ptr<ConjGradSolver> ConjGradSolver::create(const Ptr<MinProblemSolver::Function>& f, TermCriteria termcrit){
Ptr<ConjGradSolver> CG = makePtr<ConjGradSolverImpl>();
CG->setFunction(f);
CG->setTermCriteria(termcrit);
return CG;
}
}

483
modules/core/src/datastructs.cpp
View File

@@ -3528,492 +3528,9 @@ cvPrevTreeNode( CvTreeNodeIterator* treeIterator )
return prevNode;
}
namespace cv
{
// This is reimplementation of kd-trees from cvkdtree*.* by Xavier Delacour, cleaned-up and
// adopted to work with the new OpenCV data structures. It's in cxcore to be shared by
// both cv (CvFeatureTree) and ml (kNN).
// The algorithm is taken from:
// J.S. Beis and D.G. Lowe. Shape indexing using approximate nearest-neighbor search
// in highdimensional spaces. In Proc. IEEE Conf. Comp. Vision Patt. Recog.,
// pages 1000--1006, 1997. http://citeseer.ist.psu.edu/beis97shape.html
const int MAX_TREE_DEPTH = 32;
KDTree::KDTree()
{
maxDepth = -1;
normType = NORM_L2;
}
KDTree::KDTree(InputArray _points, bool _copyData)
{
maxDepth = -1;
normType = NORM_L2;
build(_points, _copyData);
}
KDTree::KDTree(InputArray _points, InputArray _labels, bool _copyData)
{
maxDepth = -1;
normType = NORM_L2;
build(_points, _labels, _copyData);
}
struct SubTree
{
SubTree() : first(0), last(0), nodeIdx(0), depth(0) {}
SubTree(int _first, int _last, int _nodeIdx, int _depth)
: first(_first), last(_last), nodeIdx(_nodeIdx), depth(_depth) {}
int first;
int last;
int nodeIdx;
int depth;
};
static float
medianPartition( size_t* ofs, int a, int b, const float* vals )
{
int k, a0 = a, b0 = b;
int middle = (a + b)/2;
while( b > a )
{
int i0 = a, i1 = (a+b)/2, i2 = b;
float v0 = vals[ofs[i0]], v1 = vals[ofs[i1]], v2 = vals[ofs[i2]];
int ip = v0 < v1 ? (v1 < v2 ? i1 : v0 < v2 ? i2 : i0) :
v0 < v2 ? i0 : (v1 < v2 ? i2 : i1);
float pivot = vals[ofs[ip]];
std::swap(ofs[ip], ofs[i2]);
for( i1 = i0, i0--; i1 <= i2; i1++ )
if( vals[ofs[i1]] <= pivot )
{
i0++;
std::swap(ofs[i0], ofs[i1]);
}
if( i0 == middle )
break;
if( i0 > middle )
b = i0 - (b == i0);
else
a = i0;
}
float pivot = vals[ofs[middle]];
int less = 0, more = 0;
for( k = a0; k < middle; k++ )
{
CV_Assert(vals[ofs[k]] <= pivot);
less += vals[ofs[k]] < pivot;
}
for( k = b0; k > middle; k-- )
{
CV_Assert(vals[ofs[k]] >= pivot);
more += vals[ofs[k]] > pivot;
}
CV_Assert(std::abs(more - less) <= 1);
return vals[ofs[middle]];
}
static void
computeSums( const Mat& points, const size_t* ofs, int a, int b, double* sums )
{
int i, j, dims = points.cols;
const float* data = points.ptr<float>(0);
for( j = 0; j < dims; j++ )
sums[j*2] = sums[j*2+1] = 0;
for( i = a; i <= b; i++ )
{
const float* row = data + ofs[i];
for( j = 0; j < dims; j++ )
{
double t = row[j], s = sums[j*2] + t, s2 = sums[j*2+1] + t*t;
sums[j*2] = s; sums[j*2+1] = s2;
}
}
}
void KDTree::build(InputArray _points, bool _copyData)
{
build(_points, noArray(), _copyData);
}
void KDTree::build(InputArray __points, InputArray __labels, bool _copyData)
{
Mat _points = __points.getMat(), _labels = __labels.getMat();
CV_Assert(_points.type() == CV_32F && !_points.empty());
std::vector<KDTree::Node>().swap(nodes);
if( !_copyData )
points = _points;
else
{
points.release();
points.create(_points.size(), _points.type());
}
int i, j, n = _points.rows, ptdims = _points.cols, top = 0;
const float* data = _points.ptr<float>(0);
float* dstdata = points.ptr<float>(0);
size_t step = _points.step1();
size_t dstep = points.step1();
int ptpos = 0;
labels.resize(n);
const int* _labels_data = 0;
if( !_labels.empty() )
{
int nlabels = _labels.checkVector(1, CV_32S, true);
CV_Assert(nlabels == n);
_labels_data = _labels.ptr<int>();
}
Mat sumstack(MAX_TREE_DEPTH*2, ptdims*2, CV_64F);
SubTree stack[MAX_TREE_DEPTH*2];
std::vector<size_t> _ptofs(n);
size_t* ptofs = &_ptofs[0];
for( i = 0; i < n; i++ )
ptofs[i] = i*step;
nodes.push_back(Node());
computeSums(points, ptofs, 0, n-1, sumstack.ptr<double>(top));
stack[top++] = SubTree(0, n-1, 0, 0);
int _maxDepth = 0;
while( --top >= 0 )
{
int first = stack[top].first, last = stack[top].last;
int depth = stack[top].depth, nidx = stack[top].nodeIdx;
int count = last - first + 1, dim = -1;
const double* sums = sumstack.ptr<double>(top);
double invCount = 1./count, maxVar = -1.;
if( count == 1 )
{
int idx0 = (int)(ptofs[first]/step);
int idx = _copyData ? ptpos++ : idx0;
nodes[nidx].idx = ~idx;
if( _copyData )
{
const float* src = data + ptofs[first];
float* dst = dstdata + idx*dstep;
for( j = 0; j < ptdims; j++ )
dst[j] = src[j];
}
labels[idx] = _labels_data ? _labels_data[idx0] : idx0;
_maxDepth = std::max(_maxDepth, depth);
continue;
}
// find the dimensionality with the biggest variance
for( j = 0; j < ptdims; j++ )
{
double m = sums[j*2]*invCount;
double varj = sums[j*2+1]*invCount - m*m;
if( maxVar < varj )
{
maxVar = varj;
dim = j;
}
}
int left = (int)nodes.size(), right = left + 1;
nodes.push_back(Node());
nodes.push_back(Node());
nodes[nidx].idx = dim;
nodes[nidx].left = left;
nodes[nidx].right = right;
nodes[nidx].boundary = medianPartition(ptofs, first, last, data + dim);
int middle = (first + last)/2;
double *lsums = (double*)sums, *rsums = lsums + ptdims*2;
computeSums(points, ptofs, middle+1, last, rsums);
for( j = 0; j < ptdims*2; j++ )
lsums[j] = sums[j] - rsums[j];
stack[top++] = SubTree(first, middle, left, depth+1);
stack[top++] = SubTree(middle+1, last, right, depth+1);
}
maxDepth = _maxDepth;
}
struct PQueueElem
{
PQueueElem() : dist(0), idx(0) {}
PQueueElem(float _dist, int _idx) : dist(_dist), idx(_idx) {}
float dist;
int idx;
};
int KDTree::findNearest(InputArray _vec, int K, int emax,
OutputArray _neighborsIdx, OutputArray _neighbors,
OutputArray _dist, OutputArray _labels) const
{
Mat vecmat = _vec.getMat();
CV_Assert( vecmat.isContinuous() && vecmat.type() == CV_32F && vecmat.total() == (size_t)points.cols );
const float* vec = vecmat.ptr<float>();
K = std::min(K, points.rows);
int ptdims = points.cols;
CV_Assert(K > 0 && (normType == NORM_L2 || normType == NORM_L1));
AutoBuffer<uchar> _buf((K+1)*(sizeof(float) + sizeof(int)));
int* idx = (int*)(uchar*)_buf;
float* dist = (float*)(idx + K + 1);
int i, j, ncount = 0, e = 0;
int qsize = 0, maxqsize = 1 << 10;
AutoBuffer<uchar> _pqueue(maxqsize*sizeof(PQueueElem));
PQueueElem* pqueue = (PQueueElem*)(uchar*)_pqueue;
emax = std::max(emax, 1);
for( e = 0; e < emax; )
{
float d, alt_d = 0.f;
int nidx;
if( e == 0 )
nidx = 0;
else
{
// take the next node from the priority queue
if( qsize == 0 )
break;
nidx = pqueue[0].idx;
alt_d = pqueue[0].dist;
if( --qsize > 0 )
{
std::swap(pqueue[0], pqueue[qsize]);
d = pqueue[0].dist;
for( i = 0;;)
{
int left = i*2 + 1, right = i*2 + 2;
if( left >= qsize )
break;
if( right < qsize && pqueue[right].dist < pqueue[left].dist )
left = right;
if( pqueue[left].dist >= d )
break;
std::swap(pqueue[i], pqueue[left]);
i = left;
}
}
if( ncount == K && alt_d > dist[ncount-1] )
continue;
}
for(;;)
{
if( nidx < 0 )
break;
const Node& n = nodes[nidx];
if( n.idx < 0 )
{
i = ~n.idx;
const float* row = points.ptr<float>(i);
if( normType == NORM_L2 )
for( j = 0, d = 0.f; j < ptdims; j++ )
{
float t = vec[j] - row[j];
d += t*t;
}
else
for( j = 0, d = 0.f; j < ptdims; j++ )
d += std::abs(vec[j] - row[j]);
dist[ncount] = d;
idx[ncount] = i;
for( i = ncount-1; i >= 0; i-- )
{
if( dist[i] <= d )
break;
std::swap(dist[i], dist[i+1]);
std::swap(idx[i], idx[i+1]);
}
ncount += ncount < K;
e++;
break;
}
int alt;
if( vec[n.idx] <= n.boundary )
{
nidx = n.left;
alt = n.right;
}
else
{
nidx = n.right;
alt = n.left;
}
d = vec[n.idx] - n.boundary;
if( normType == NORM_L2 )
d = d*d + alt_d;
else
d = std::abs(d) + alt_d;
// subtree prunning
if( ncount == K && d > dist[ncount-1] )
continue;
// add alternative subtree to the priority queue
pqueue[qsize] = PQueueElem(d, alt);
for( i = qsize; i > 0; )
{
int parent = (i-1)/2;
if( parent < 0 || pqueue[parent].dist <= d )
break;
std::swap(pqueue[i], pqueue[parent]);
i = parent;
}
qsize += qsize+1 < maxqsize;
}
}
K = std::min(K, ncount);
if( _neighborsIdx.needed() )
{
_neighborsIdx.create(K, 1, CV_32S, -1, true);
Mat nidx = _neighborsIdx.getMat();
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
}
if( _dist.needed() )
sqrt(Mat(K, 1, CV_32F, dist), _dist);
if( _neighbors.needed() || _labels.needed() )
getPoints(Mat(K, 1, CV_32S, idx), _neighbors, _labels);
return K;
}
void KDTree::findOrthoRange(InputArray _lowerBound,
InputArray _upperBound,
OutputArray _neighborsIdx,
OutputArray _neighbors,
OutputArray _labels ) const
{
int ptdims = points.cols;
Mat lowerBound = _lowerBound.getMat(), upperBound = _upperBound.getMat();
CV_Assert( lowerBound.size == upperBound.size &&
lowerBound.isContinuous() &&
upperBound.isContinuous() &&
lowerBound.type() == upperBound.type() &&
lowerBound.type() == CV_32F &&
lowerBound.total() == (size_t)ptdims );
const float* L = lowerBound.ptr<float>();
const float* R = upperBound.ptr<float>();
std::vector<int> idx;
AutoBuffer<int> _stack(MAX_TREE_DEPTH*2 + 1);
int* stack = _stack;
int top = 0;
stack[top++] = 0;
while( --top >= 0 )
{
int nidx = stack[top];
if( nidx < 0 )
break;
const Node& n = nodes[nidx];
if( n.idx < 0 )
{
int j, i = ~n.idx;
const float* row = points.ptr<float>(i);
for( j = 0; j < ptdims; j++ )
if( row[j] < L[j] || row[j] >= R[j] )
break;
if( j == ptdims )
idx.push_back(i);
continue;
}
if( L[n.idx] <= n.boundary )
stack[top++] = n.left;
if( R[n.idx] > n.boundary )
stack[top++] = n.right;
}
if( _neighborsIdx.needed() )
{
_neighborsIdx.create((int)idx.size(), 1, CV_32S, -1, true);
Mat nidx = _neighborsIdx.getMat();
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
}
getPoints( idx, _neighbors, _labels );
}
void KDTree::getPoints(InputArray _idx, OutputArray _pts, OutputArray _labels) const
{
Mat idxmat = _idx.getMat(), pts, labelsmat;
CV_Assert( idxmat.isContinuous() && idxmat.type() == CV_32S &&
(idxmat.cols == 1 || idxmat.rows == 1) );
const int* idx = idxmat.ptr<int>();
int* dstlabels = 0;
int ptdims = points.cols;
int i, nidx = (int)idxmat.total();
if( nidx == 0 )
{
_pts.release();
_labels.release();
return;
}
if( _pts.needed() )
{
_pts.create( nidx, ptdims, points.type());
pts = _pts.getMat();
}
if(_labels.needed())
{
_labels.create(nidx, 1, CV_32S, -1, true);
labelsmat = _labels.getMat();
CV_Assert( labelsmat.isContinuous() );
dstlabels = labelsmat.ptr<int>();
}
const int* srclabels = !labels.empty() ? &labels[0] : 0;
for( i = 0; i < nidx; i++ )
{
int k = idx[i];
CV_Assert( (unsigned)k < (unsigned)points.rows );
const float* src = points.ptr<float>(k);
if( !pts.empty() )
std::copy(src, src + ptdims, pts.ptr<float>(i));
if( dstlabels )
dstlabels[i] = srclabels ? srclabels[k] : k;
}
}
const float* KDTree::getPoint(int ptidx, int* label) const
{
CV_Assert( (unsigned)ptidx < (unsigned)points.rows);
if(label)
*label = labels[ptidx];
return points.ptr<float>(ptidx);
}
int KDTree::dims() const
{
return !points.empty() ? points.cols : 0;
}
////////////////////////////////////////////////////////////////////////////////
schar* seqPush( CvSeq* seq, const void* element )

409
modules/core/src/downhill_simplex.cpp Normal file
View File

@@ -0,0 +1,409 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#define dprintf(x)
#define print_matrix(x)
/*
****Error Message********************************************************************************************************************
Downhill Simplex method in OpenCV dev 3.0.0 getting this error:
OpenCV Error: Assertion failed (dims <= 2 && data && (unsigned)i0 < (unsigned)(s ize.p[0] * size.p[1])
&& elemSize() == (((((DataType<_Tp>::type) & ((512 - 1) << 3)) >> 3) + 1) << ((((sizeof(size_t)/4+1)16384|0x3a50)
>> ((DataType<_Tp>::typ e) & ((1 << 3) - 1))2) & 3))) in cv::Mat::at,
file C:\builds\master_PackSlave-w in32-vc12-shared\opencv\modules\core\include\opencv2/core/mat.inl.hpp, line 893
****Problem and Possible Fix*********************************************************************************************************
DownhillSolverImpl::innerDownhillSimplex something looks broken here:
Mat_<double> coord_sum(1,ndim,0.0),buf(1,ndim,0.0),y(1,ndim,0.0);
nfunk = 0;
for(i=0;i<ndim+1;++i)
{
y(i) = f->calc(p[i]);
}
y has only ndim elements, while the loop goes over ndim+1
Edited the following for possible fix:
Replaced y(1,ndim,0.0) ------> y(1,ndim+1,0.0)
***********************************************************************************************************************************
The code below was used in tesing the source code.
Created by @SareeAlnaghy
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <opencv2\optim\optim.hpp>
using namespace std;
using namespace cv;
void test(Ptr<optim::DownhillSolver> MinProblemSolver, Ptr<optim::MinProblemSolver::Function> ptr_F, Mat &P, Mat &step)
{
try{
MinProblemSolver->setFunction(ptr_F);
MinProblemSolver->setInitStep(step);
double res = MinProblemSolver->minimize(P);
cout << "res " << res << endl;
}
catch (exception e)
{
cerr << "Error:: " << e.what() << endl;
}
}
int main()
{
class DistanceToLines :public optim::MinProblemSolver::Function {
public:
double calc(const double* x)const{
return x[0] * x[0] + x[1] * x[1];
}
};
Mat P = (Mat_<double>(1, 2) << 1.0, 1.0);
Mat step = (Mat_<double>(2, 1) << -0.5, 0.5);
Ptr<optim::MinProblemSolver::Function> ptr_F(new DistanceToLines());
Ptr<optim::DownhillSolver> MinProblemSolver = optim::createDownhillSolver();
test(MinProblemSolver, ptr_F, P, step);
system("pause");
return 0;
}
****Suggesttion for imporving Simplex implentation***************************************************************************************
Currently the downhilll simplex method outputs the function value that is minimized. It should also return the coordinate points where the
function is minimized. This is very useful in many applications such as using back projection methods to find a point of intersection of
multiple lines in three dimensions as not all lines intersect in three dimensions.
*/
namespace cv
{
class DownhillSolverImpl : public DownhillSolver
{
public:
void getInitStep(OutputArray step) const;
void setInitStep(InputArray step);
Ptr<Function> getFunction() const;
void setFunction(const Ptr<Function>& f);
TermCriteria getTermCriteria() const;
DownhillSolverImpl();
void setTermCriteria(const TermCriteria& termcrit);
double minimize(InputOutputArray x);
protected:
Ptr<MinProblemSolver::Function> _Function;
TermCriteria _termcrit;
Mat _step;
Mat_<double> buf_x;
private:
inline void createInitialSimplex(Mat_<double>& simplex,Mat& step);
inline double innerDownhillSimplex(cv::Mat_<double>& p,double MinRange,double MinError,int& nfunk,
const Ptr<MinProblemSolver::Function>& f,int nmax);
inline double tryNewPoint(Mat_<double>& p,Mat_<double>& y,Mat_<double>& coord_sum,const Ptr<MinProblemSolver::Function>& f,int ihi,
double fac,Mat_<double>& ptry);
};
double DownhillSolverImpl::tryNewPoint(
Mat_<double>& p,
Mat_<double>& y,
Mat_<double>& coord_sum,
const Ptr<MinProblemSolver::Function>& f,
int ihi,
double fac,
Mat_<double>& ptry
)
{
int ndim=p.cols;
int j;
double fac1,fac2,ytry;
fac1=(1.0-fac)/ndim;
fac2=fac1-fac;
for (j=0;j<ndim;j++)
{
ptry(j)=coord_sum(j)*fac1-p(ihi,j)*fac2;
}
ytry=f->calc(ptry.ptr<double>());
if (ytry < y(ihi))
{
y(ihi)=ytry;
for (j=0;j<ndim;j++)
{
coord_sum(j) += ptry(j)-p(ihi,j);
p(ihi,j)=ptry(j);
}
}
return ytry;
}
/*
Performs the actual minimization of MinProblemSolver::Function f (after the initialization was done)
The matrix p[ndim+1][1..ndim] represents ndim+1 vertices that
form a simplex - each row is an ndim vector.
On output, nfunk gives the number of function evaluations taken.
*/
double DownhillSolverImpl::innerDownhillSimplex(
cv::Mat_<double>& p,
double MinRange,
double MinError,
int& nfunk,
const Ptr<MinProblemSolver::Function>& f,
int nmax
)
{
int ndim=p.cols;
double res;
int i,ihi,ilo,inhi,j,mpts=ndim+1;
double error, range,ysave,ytry;
Mat_<double> coord_sum(1,ndim,0.0),buf(1,ndim,0.0),y(1,ndim+1,0.0);
nfunk = 0;
for(i=0;i<ndim+1;++i)
{
y(i) = f->calc(p[i]);
}
nfunk = ndim+1;
reduce(p,coord_sum,0,CV_REDUCE_SUM);
for (;;)
{
ilo=0;
/* find highest (worst), next-to-worst, and lowest
(best) points by going through all of them. */
ihi = y(0)>y(1) ? (inhi=1,0) : (inhi=0,1);
for (i=0;i<mpts;i++)
{
if (y(i) <= y(ilo))
ilo=i;
if (y(i) > y(ihi))
{
inhi=ihi;
ihi=i;
}
else if (y(i) > y(inhi) && i != ihi)
inhi=i;
}
/* check stop criterion */
error=fabs(y(ihi)-y(ilo));
range=0;
for(i=0;i<ndim;++i)
{
double min = p(0,i);
double max = p(0,i);
double d;
for(j=1;j<=ndim;++j)
{
if( min > p(j,i) ) min = p(j,i);
if( max < p(j,i) ) max = p(j,i);
}
d = fabs(max-min);
if(range < d) range = d;
}
if(range <= MinRange || error <= MinError)
{ /* Put best point and value in first slot. */
std::swap(y(0),y(ilo));
for (i=0;i<ndim;i++)
{
std::swap(p(0,i),p(ilo,i));
}
break;
}
if (nfunk >= nmax){
dprintf(("nmax exceeded\n"));
return y(ilo);
}
nfunk += 2;
/*Begin a new iteration. First, reflect the worst point about the centroid of others */
ytry = tryNewPoint(p,y,coord_sum,f,ihi,-1.0,buf);
if (ytry <= y(ilo))
{ /*If that's better than the best point, go twice as far in that direction*/
ytry = tryNewPoint(p,y,coord_sum,f,ihi,2.0,buf);
}
else if (ytry >= y(inhi))
{ /* The new point is worse than the second-highest, but better
than the worst so do not go so far in that direction */
ysave = y(ihi);
ytry = tryNewPoint(p,y,coord_sum,f,ihi,0.5,buf);
if (ytry >= ysave)
{ /* Can't seem to improve things. Contract the simplex to good point
in hope to find a simplex landscape. */
for (i=0;i<mpts;i++)
{
if (i != ilo)
{
for (j=0;j<ndim;j++)
{
p(i,j) = coord_sum(j) = 0.5*(p(i,j)+p(ilo,j));
}
y(i)=f->calc(coord_sum.ptr<double>());
}
}
nfunk += ndim;
reduce(p,coord_sum,0,CV_REDUCE_SUM);
}
} else --(nfunk); /* correct nfunk */
dprintf(("this is simplex on iteration %d\n",nfunk));
print_matrix(p);
} /* go to next iteration. */
res = y(0);
return res;
}
void DownhillSolverImpl::createInitialSimplex(Mat_<double>& simplex,Mat& step){
for(int i=1;i<=step.cols;++i)
{
simplex.row(0).copyTo(simplex.row(i));
simplex(i,i-1)+= 0.5*step.at<double>(0,i-1);
}
simplex.row(0) -= 0.5*step;
dprintf(("this is simplex\n"));
print_matrix(simplex);
}
double DownhillSolverImpl::minimize(InputOutputArray x){
dprintf(("hi from minimize\n"));
CV_Assert(_Function.empty()==false);
dprintf(("termcrit:\n\ttype: %d\n\tmaxCount: %d\n\tEPS: %g\n",_termcrit.type,_termcrit.maxCount,_termcrit.epsilon));
dprintf(("step\n"));
print_matrix(_step);
Mat x_mat=x.getMat();
CV_Assert(MIN(x_mat.rows,x_mat.cols)==1);
CV_Assert(MAX(x_mat.rows,x_mat.cols)==_step.cols);
CV_Assert(x_mat.type()==CV_64FC1);
Mat_<double> proxy_x;
if(x_mat.rows>1){
buf_x.create(1,_step.cols);
Mat_<double> proxy(_step.cols,1,buf_x.ptr<double>());
x_mat.copyTo(proxy);
proxy_x=buf_x;
}else{
proxy_x=x_mat;
}
int count=0;
int ndim=_step.cols;
Mat_<double> simplex=Mat_<double>(ndim+1,ndim,0.0);
simplex.row(0).copyTo(proxy_x);
createInitialSimplex(simplex,_step);
double res = innerDownhillSimplex(
simplex,_termcrit.epsilon, _termcrit.epsilon, count,_Function,_termcrit.maxCount);
simplex.row(0).copyTo(proxy_x);
dprintf(("%d iterations done\n",count));
if(x_mat.rows>1){
Mat(x_mat.rows, 1, CV_64F, proxy_x.ptr<double>()).copyTo(x);
}
return res;
}
DownhillSolverImpl::DownhillSolverImpl(){
_Function=Ptr<Function>();
_step=Mat_<double>();
}
Ptr<MinProblemSolver::Function> DownhillSolverImpl::getFunction()const{
return _Function;
}
void DownhillSolverImpl::setFunction(const Ptr<Function>& f){
_Function=f;
}
TermCriteria DownhillSolverImpl::getTermCriteria()const{
return _termcrit;
}
void DownhillSolverImpl::setTermCriteria(const TermCriteria& termcrit){
CV_Assert(termcrit.type==(TermCriteria::MAX_ITER+TermCriteria::EPS) && termcrit.epsilon>0 && termcrit.maxCount>0);
_termcrit=termcrit;
}
// both minRange & minError are specified by termcrit.epsilon; In addition, user may specify the number of iterations that the algorithm does.
Ptr<DownhillSolver> DownhillSolver::create(const Ptr<MinProblemSolver::Function>& f, InputArray initStep, TermCriteria termcrit){
Ptr<DownhillSolver> DS = makePtr<DownhillSolverImpl>();
DS->setFunction(f);
DS->setInitStep(initStep);
DS->setTermCriteria(termcrit);
return DS;
}
void DownhillSolverImpl::getInitStep(OutputArray step)const{
_step.copyTo(step);
}
void DownhillSolverImpl::setInitStep(InputArray step){
//set dimensionality and make a deep copy of step
Mat m=step.getMat();
dprintf(("m.cols=%d\nm.rows=%d\n",m.cols,m.rows));
CV_Assert(MIN(m.cols,m.rows)==1 && m.type()==CV_64FC1);
if(m.rows==1){
m.copyTo(_step);
}else{
transpose(m,_step);
}
}
}

2552
modules/core/src/drawing.cpp
View File

File diff suppressed because it is too large Load Diff

531
modules/core/src/kdtree.cpp Normal file
View File

@@ -0,0 +1,531 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
namespace cv
{
// This is reimplementation of kd-trees from cvkdtree*.* by Xavier Delacour, cleaned-up and
// adopted to work with the new OpenCV data structures. It's in cxcore to be shared by
// both cv (CvFeatureTree) and ml (kNN).
// The algorithm is taken from:
// J.S. Beis and D.G. Lowe. Shape indexing using approximate nearest-neighbor search
// in highdimensional spaces. In Proc. IEEE Conf. Comp. Vision Patt. Recog.,
// pages 1000--1006, 1997. http://citeseer.ist.psu.edu/beis97shape.html
const int MAX_TREE_DEPTH = 32;
KDTree::KDTree()
{
maxDepth = -1;
normType = NORM_L2;
}
KDTree::KDTree(InputArray _points, bool _copyData)
{
maxDepth = -1;
normType = NORM_L2;
build(_points, _copyData);
}
KDTree::KDTree(InputArray _points, InputArray _labels, bool _copyData)
{
maxDepth = -1;
normType = NORM_L2;
build(_points, _labels, _copyData);
}
struct SubTree
{
SubTree() : first(0), last(0), nodeIdx(0), depth(0) {}
SubTree(int _first, int _last, int _nodeIdx, int _depth)
: first(_first), last(_last), nodeIdx(_nodeIdx), depth(_depth) {}
int first;
int last;
int nodeIdx;
int depth;
};
static float
medianPartition( size_t* ofs, int a, int b, const float* vals )
{
int k, a0 = a, b0 = b;
int middle = (a + b)/2;
while( b > a )
{
int i0 = a, i1 = (a+b)/2, i2 = b;
float v0 = vals[ofs[i0]], v1 = vals[ofs[i1]], v2 = vals[ofs[i2]];
int ip = v0 < v1 ? (v1 < v2 ? i1 : v0 < v2 ? i2 : i0) :
v0 < v2 ? i0 : (v1 < v2 ? i2 : i1);
float pivot = vals[ofs[ip]];
std::swap(ofs[ip], ofs[i2]);
for( i1 = i0, i0--; i1 <= i2; i1++ )
if( vals[ofs[i1]] <= pivot )
{
i0++;
std::swap(ofs[i0], ofs[i1]);
}
if( i0 == middle )
break;
if( i0 > middle )
b = i0 - (b == i0);
else
a = i0;
}
float pivot = vals[ofs[middle]];
int less = 0, more = 0;
for( k = a0; k < middle; k++ )
{
CV_Assert(vals[ofs[k]] <= pivot);
less += vals[ofs[k]] < pivot;
}
for( k = b0; k > middle; k-- )
{
CV_Assert(vals[ofs[k]] >= pivot);
more += vals[ofs[k]] > pivot;
}
CV_Assert(std::abs(more - less) <= 1);
return vals[ofs[middle]];
}
static void
computeSums( const Mat& points, const size_t* ofs, int a, int b, double* sums )
{
int i, j, dims = points.cols;
const float* data = points.ptr<float>(0);
for( j = 0; j < dims; j++ )
sums[j*2] = sums[j*2+1] = 0;
for( i = a; i <= b; i++ )
{
const float* row = data + ofs[i];
for( j = 0; j < dims; j++ )
{
double t = row[j], s = sums[j*2] + t, s2 = sums[j*2+1] + t*t;
sums[j*2] = s; sums[j*2+1] = s2;
}
}
}
void KDTree::build(InputArray _points, bool _copyData)
{
build(_points, noArray(), _copyData);
}
void KDTree::build(InputArray __points, InputArray __labels, bool _copyData)
{
Mat _points = __points.getMat(), _labels = __labels.getMat();
CV_Assert(_points.type() == CV_32F && !_points.empty());
std::vector<KDTree::Node>().swap(nodes);
if( !_copyData )
points = _points;
else
{
points.release();
points.create(_points.size(), _points.type());
}
int i, j, n = _points.rows, ptdims = _points.cols, top = 0;
const float* data = _points.ptr<float>(0);
float* dstdata = points.ptr<float>(0);
size_t step = _points.step1();
size_t dstep = points.step1();
int ptpos = 0;
labels.resize(n);
const int* _labels_data = 0;
if( !_labels.empty() )
{
int nlabels = _labels.checkVector(1, CV_32S, true);
CV_Assert(nlabels == n);
_labels_data = _labels.ptr<int>();
}
Mat sumstack(MAX_TREE_DEPTH*2, ptdims*2, CV_64F);
SubTree stack[MAX_TREE_DEPTH*2];
std::vector<size_t> _ptofs(n);
size_t* ptofs = &_ptofs[0];
for( i = 0; i < n; i++ )
ptofs[i] = i*step;
nodes.push_back(Node());
computeSums(points, ptofs, 0, n-1, sumstack.ptr<double>(top));
stack[top++] = SubTree(0, n-1, 0, 0);
int _maxDepth = 0;
while( --top >= 0 )
{
int first = stack[top].first, last = stack[top].last;
int depth = stack[top].depth, nidx = stack[top].nodeIdx;
int count = last - first + 1, dim = -1;
const double* sums = sumstack.ptr<double>(top);
double invCount = 1./count, maxVar = -1.;
if( count == 1 )
{
int idx0 = (int)(ptofs[first]/step);
int idx = _copyData ? ptpos++ : idx0;
nodes[nidx].idx = ~idx;
if( _copyData )
{
const float* src = data + ptofs[first];
float* dst = dstdata + idx*dstep;
for( j = 0; j < ptdims; j++ )
dst[j] = src[j];
}
labels[idx] = _labels_data ? _labels_data[idx0] : idx0;
_maxDepth = std::max(_maxDepth, depth);
continue;
}
// find the dimensionality with the biggest variance
for( j = 0; j < ptdims; j++ )
{
double m = sums[j*2]*invCount;
double varj = sums[j*2+1]*invCount - m*m;
if( maxVar < varj )
{
maxVar = varj;
dim = j;
}
}
int left = (int)nodes.size(), right = left + 1;
nodes.push_back(Node());
nodes.push_back(Node());
nodes[nidx].idx = dim;
nodes[nidx].left = left;
nodes[nidx].right = right;
nodes[nidx].boundary = medianPartition(ptofs, first, last, data + dim);
int middle = (first + last)/2;
double *lsums = (double*)sums, *rsums = lsums + ptdims*2;
computeSums(points, ptofs, middle+1, last, rsums);
for( j = 0; j < ptdims*2; j++ )
lsums[j] = sums[j] - rsums[j];
stack[top++] = SubTree(first, middle, left, depth+1);
stack[top++] = SubTree(middle+1, last, right, depth+1);
}
maxDepth = _maxDepth;
}
struct PQueueElem
{
PQueueElem() : dist(0), idx(0) {}
PQueueElem(float _dist, int _idx) : dist(_dist), idx(_idx) {}
float dist;
int idx;
};
int KDTree::findNearest(InputArray _vec, int K, int emax,
OutputArray _neighborsIdx, OutputArray _neighbors,
OutputArray _dist, OutputArray _labels) const
{
Mat vecmat = _vec.getMat();
CV_Assert( vecmat.isContinuous() && vecmat.type() == CV_32F && vecmat.total() == (size_t)points.cols );
const float* vec = vecmat.ptr<float>();
K = std::min(K, points.rows);
int ptdims = points.cols;
CV_Assert(K > 0 && (normType == NORM_L2 || normType == NORM_L1));
AutoBuffer<uchar> _buf((K+1)*(sizeof(float) + sizeof(int)));
int* idx = (int*)(uchar*)_buf;
float* dist = (float*)(idx + K + 1);
int i, j, ncount = 0, e = 0;
int qsize = 0, maxqsize = 1 << 10;
AutoBuffer<uchar> _pqueue(maxqsize*sizeof(PQueueElem));
PQueueElem* pqueue = (PQueueElem*)(uchar*)_pqueue;
emax = std::max(emax, 1);
for( e = 0; e < emax; )
{
float d, alt_d = 0.f;
int nidx;
if( e == 0 )
nidx = 0;
else
{
// take the next node from the priority queue
if( qsize == 0 )
break;
nidx = pqueue[0].idx;
alt_d = pqueue[0].dist;
if( --qsize > 0 )
{
std::swap(pqueue[0], pqueue[qsize]);
d = pqueue[0].dist;
for( i = 0;;)
{
int left = i*2 + 1, right = i*2 + 2;
if( left >= qsize )
break;
if( right < qsize && pqueue[right].dist < pqueue[left].dist )
left = right;
if( pqueue[left].dist >= d )
break;
std::swap(pqueue[i], pqueue[left]);
i = left;
}
}
if( ncount == K && alt_d > dist[ncount-1] )
continue;
}
for(;;)
{
if( nidx < 0 )
break;
const Node& n = nodes[nidx];
if( n.idx < 0 )
{
i = ~n.idx;
const float* row = points.ptr<float>(i);
if( normType == NORM_L2 )
for( j = 0, d = 0.f; j < ptdims; j++ )
{
float t = vec[j] - row[j];
d += t*t;
}
else
for( j = 0, d = 0.f; j < ptdims; j++ )
d += std::abs(vec[j] - row[j]);
dist[ncount] = d;
idx[ncount] = i;
for( i = ncount-1; i >= 0; i-- )
{
if( dist[i] <= d )
break;
std::swap(dist[i], dist[i+1]);
std::swap(idx[i], idx[i+1]);
}
ncount += ncount < K;
e++;
break;
}
int alt;
if( vec[n.idx] <= n.boundary )
{
nidx = n.left;
alt = n.right;
}
else
{
nidx = n.right;
alt = n.left;
}
d = vec[n.idx] - n.boundary;
if( normType == NORM_L2 )
d = d*d + alt_d;
else
d = std::abs(d) + alt_d;
// subtree prunning
if( ncount == K && d > dist[ncount-1] )
continue;
// add alternative subtree to the priority queue
pqueue[qsize] = PQueueElem(d, alt);
for( i = qsize; i > 0; )
{
int parent = (i-1)/2;
if( parent < 0 || pqueue[parent].dist <= d )
break;
std::swap(pqueue[i], pqueue[parent]);
i = parent;
}
qsize += qsize+1 < maxqsize;
}
}
K = std::min(K, ncount);
if( _neighborsIdx.needed() )
{
_neighborsIdx.create(K, 1, CV_32S, -1, true);
Mat nidx = _neighborsIdx.getMat();
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
}
if( _dist.needed() )
sqrt(Mat(K, 1, CV_32F, dist), _dist);
if( _neighbors.needed() || _labels.needed() )
getPoints(Mat(K, 1, CV_32S, idx), _neighbors, _labels);
return K;
}
void KDTree::findOrthoRange(InputArray _lowerBound,
InputArray _upperBound,
OutputArray _neighborsIdx,
OutputArray _neighbors,
OutputArray _labels ) const
{
int ptdims = points.cols;
Mat lowerBound = _lowerBound.getMat(), upperBound = _upperBound.getMat();
CV_Assert( lowerBound.size == upperBound.size &&
lowerBound.isContinuous() &&
upperBound.isContinuous() &&
lowerBound.type() == upperBound.type() &&
lowerBound.type() == CV_32F &&
lowerBound.total() == (size_t)ptdims );
const float* L = lowerBound.ptr<float>();
const float* R = upperBound.ptr<float>();
std::vector<int> idx;
AutoBuffer<int> _stack(MAX_TREE_DEPTH*2 + 1);
int* stack = _stack;
int top = 0;
stack[top++] = 0;
while( --top >= 0 )
{
int nidx = stack[top];
if( nidx < 0 )
break;
const Node& n = nodes[nidx];
if( n.idx < 0 )
{
int j, i = ~n.idx;
const float* row = points.ptr<float>(i);
for( j = 0; j < ptdims; j++ )
if( row[j] < L[j] || row[j] >= R[j] )
break;
if( j == ptdims )
idx.push_back(i);
continue;
}
if( L[n.idx] <= n.boundary )
stack[top++] = n.left;
if( R[n.idx] > n.boundary )
stack[top++] = n.right;
}
if( _neighborsIdx.needed() )
{
_neighborsIdx.create((int)idx.size(), 1, CV_32S, -1, true);
Mat nidx = _neighborsIdx.getMat();
Mat(nidx.size(), CV_32S, &idx[0]).copyTo(nidx);
}
getPoints( idx, _neighbors, _labels );
}
void KDTree::getPoints(InputArray _idx, OutputArray _pts, OutputArray _labels) const
{
Mat idxmat = _idx.getMat(), pts, labelsmat;
CV_Assert( idxmat.isContinuous() && idxmat.type() == CV_32S &&
(idxmat.cols == 1 || idxmat.rows == 1) );
const int* idx = idxmat.ptr<int>();
int* dstlabels = 0;
int ptdims = points.cols;
int i, nidx = (int)idxmat.total();
if( nidx == 0 )
{
_pts.release();
_labels.release();
return;
}
if( _pts.needed() )
{
_pts.create( nidx, ptdims, points.type());
pts = _pts.getMat();
}
if(_labels.needed())
{
_labels.create(nidx, 1, CV_32S, -1, true);
labelsmat = _labels.getMat();
CV_Assert( labelsmat.isContinuous() );
dstlabels = labelsmat.ptr<int>();
}
const int* srclabels = !labels.empty() ? &labels[0] : 0;
for( i = 0; i < nidx; i++ )
{
int k = idx[i];
CV_Assert( (unsigned)k < (unsigned)points.rows );
const float* src = points.ptr<float>(k);
if( !pts.empty() )
std::copy(src, src + ptdims, pts.ptr<float>(i));
if( dstlabels )
dstlabels[i] = srclabels ? srclabels[k] : k;
}
}
const float* KDTree::getPoint(int ptidx, int* label) const
{
CV_Assert( (unsigned)ptidx < (unsigned)points.rows);
if(label)
*label = labels[ptidx];
return points.ptr<float>(ptidx);
}
int KDTree::dims() const
{
return !points.empty() ? points.cols : 0;
}
}

457
modules/core/src/kmeans.cpp Normal file
View File

@@ -0,0 +1,457 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
////////////////////////////////////////// kmeans ////////////////////////////////////////////
namespace cv
{
static void generateRandomCenter(const std::vector<Vec2f>& box, float* center, RNG& rng)
{
size_t j, dims = box.size();
float margin = 1.f/dims;
for( j = 0; j < dims; j++ )
center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
}
class KMeansPPDistanceComputer : public ParallelLoopBody
{
public:
KMeansPPDistanceComputer( float *_tdist2,
const float *_data,
const float *_dist,
int _dims,
size_t _step,
size_t _stepci )
: tdist2(_tdist2),
data(_data),
dist(_dist),
dims(_dims),
step(_step),
stepci(_stepci) { }
void operator()( const cv::Range& range ) const
{
const int begin = range.start;
const int end = range.end;
for ( int i = begin; i<end; i++ )
{
tdist2[i] = std::min(normL2Sqr_(data + step*i, data + stepci, dims), dist[i]);
}
}
private:
KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // to quiet MSVC
float *tdist2;
const float *data;
const float *dist;
const int dims;
const size_t step;
const size_t stepci;
};
/*
k-means center initialization using the following algorithm:
Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding
*/
static void generateCentersPP(const Mat& _data, Mat& _out_centers,
int K, RNG& rng, int trials)
{
int i, j, k, dims = _data.cols, N = _data.rows;
const float* data = _data.ptr<float>(0);
size_t step = _data.step/sizeof(data[0]);
std::vector<int> _centers(K);
int* centers = &_centers[0];
std::vector<float> _dist(N*3);
float* dist = &_dist[0], *tdist = dist + N, *tdist2 = tdist + N;
double sum0 = 0;
centers[0] = (unsigned)rng % N;
for( i = 0; i < N; i++ )
{
dist[i] = normL2Sqr_(data + step*i, data + step*centers[0], dims);
sum0 += dist[i];
}
for( k = 1; k < K; k++ )
{
double bestSum = DBL_MAX;
int bestCenter = -1;
for( j = 0; j < trials; j++ )
{
double p = (double)rng*sum0, s = 0;
for( i = 0; i < N-1; i++ )
if( (p -= dist[i]) <= 0 )
break;
int ci = i;
parallel_for_(Range(0, N),
KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
for( i = 0; i < N; i++ )
{
s += tdist2[i];
}
if( s < bestSum )
{
bestSum = s;
bestCenter = ci;
std::swap(tdist, tdist2);
}
}
centers[k] = bestCenter;
sum0 = bestSum;
std::swap(dist, tdist);
}
for( k = 0; k < K; k++ )
{
const float* src = data + step*centers[k];
float* dst = _out_centers.ptr<float>(k);
for( j = 0; j < dims; j++ )
dst[j] = src[j];
}
}
class KMeansDistanceComputer : public ParallelLoopBody
{
public:
KMeansDistanceComputer( double *_distances,
int *_labels,
const Mat& _data,
const Mat& _centers )
: distances(_distances),
labels(_labels),
data(_data),
centers(_centers)
{
}
void operator()( const Range& range ) const
{
const int begin = range.start;
const int end = range.end;
const int K = centers.rows;
const int dims = centers.cols;
const float *sample;
for( int i = begin; i<end; ++i)
{
sample = data.ptr<float>(i);
int k_best = 0;
double min_dist = DBL_MAX;
for( int k = 0; k < K; k++ )
{
const float* center = centers.ptr<float>(k);
const double dist = normL2Sqr_(sample, center, dims);
if( min_dist > dist )
{
min_dist = dist;
k_best = k;
}
}
distances[i] = min_dist;
labels[i] = k_best;
}
}
private:
KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // to quiet MSVC
double *distances;
int *labels;
const Mat& data;
const Mat& centers;
};
}
double cv::kmeans( InputArray _data, int K,
InputOutputArray _bestLabels,
TermCriteria criteria, int attempts,
int flags, OutputArray _centers )
{
const int SPP_TRIALS = 3;
Mat data0 = _data.getMat();
bool isrow = data0.rows == 1 && data0.channels() > 1;
int N = !isrow ? data0.rows : data0.cols;
int dims = (!isrow ? data0.cols : 1)*data0.channels();
int type = data0.depth();
attempts = std::max(attempts, 1);
CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
CV_Assert( N >= K );
Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));
_bestLabels.create(N, 1, CV_32S, -1, true);
Mat _labels, best_labels = _bestLabels.getMat();
if( flags & CV_KMEANS_USE_INITIAL_LABELS )
{
CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
best_labels.cols*best_labels.rows == N &&
best_labels.type() == CV_32S &&
best_labels.isContinuous());
best_labels.copyTo(_labels);
}
else
{
if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
best_labels.cols*best_labels.rows == N &&
best_labels.type() == CV_32S &&
best_labels.isContinuous()))
best_labels.create(N, 1, CV_32S);
_labels.create(best_labels.size(), best_labels.type());
}
int* labels = _labels.ptr<int>();
Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
std::vector<int> counters(K);
std::vector<Vec2f> _box(dims);
Vec2f* box = &_box[0];
double best_compactness = DBL_MAX, compactness = 0;
RNG& rng = theRNG();
int a, iter, i, j, k;
if( criteria.type & TermCriteria::EPS )
criteria.epsilon = std::max(criteria.epsilon, 0.);
else
criteria.epsilon = FLT_EPSILON;
criteria.epsilon *= criteria.epsilon;
if( criteria.type & TermCriteria::COUNT )
criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
else
criteria.maxCount = 100;
if( K == 1 )
{
attempts = 1;
criteria.maxCount = 2;
}
const float* sample = data.ptr<float>(0);
for( j = 0; j < dims; j++ )
box[j] = Vec2f(sample[j], sample[j]);
for( i = 1; i < N; i++ )
{
sample = data.ptr<float>(i);
for( j = 0; j < dims; j++ )
{
float v = sample[j];
box[j][0] = std::min(box[j][0], v);
box[j][1] = std::max(box[j][1], v);
}
}
for( a = 0; a < attempts; a++ )
{
double max_center_shift = DBL_MAX;
for( iter = 0;; )
{
swap(centers, old_centers);
if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
{
if( flags & KMEANS_PP_CENTERS )
generateCentersPP(data, centers, K, rng, SPP_TRIALS);
else
{
for( k = 0; k < K; k++ )
generateRandomCenter(_box, centers.ptr<float>(k), rng);
}
}
else
{
if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
{
for( i = 0; i < N; i++ )
CV_Assert( (unsigned)labels[i] < (unsigned)K );
}
// compute centers
centers = Scalar(0);
for( k = 0; k < K; k++ )
counters[k] = 0;
for( i = 0; i < N; i++ )
{
sample = data.ptr<float>(i);
k = labels[i];
float* center = centers.ptr<float>(k);
j=0;
#if CV_ENABLE_UNROLLED
for(; j <= dims - 4; j += 4 )
{
float t0 = center[j] + sample[j];
float t1 = center[j+1] + sample[j+1];
center[j] = t0;
center[j+1] = t1;
t0 = center[j+2] + sample[j+2];
t1 = center[j+3] + sample[j+3];
center[j+2] = t0;
center[j+3] = t1;
}
#endif
for( ; j < dims; j++ )
center[j] += sample[j];
counters[k]++;
}
if( iter > 0 )
max_center_shift = 0;
for( k = 0; k < K; k++ )
{
if( counters[k] != 0 )
continue;
// if some cluster appeared to be empty then:
// 1. find the biggest cluster
// 2. find the farthest from the center point in the biggest cluster
// 3. exclude the farthest point from the biggest cluster and form a new 1-point cluster.
int max_k = 0;
for( int k1 = 1; k1 < K; k1++ )
{
if( counters[max_k] < counters[k1] )
max_k = k1;
}
double max_dist = 0;
int farthest_i = -1;
float* new_center = centers.ptr<float>(k);
float* old_center = centers.ptr<float>(max_k);
float* _old_center = temp.ptr<float>(); // normalized
float scale = 1.f/counters[max_k];
for( j = 0; j < dims; j++ )
_old_center[j] = old_center[j]*scale;
for( i = 0; i < N; i++ )
{
if( labels[i] != max_k )
continue;
sample = data.ptr<float>(i);
double dist = normL2Sqr_(sample, _old_center, dims);
if( max_dist <= dist )
{
max_dist = dist;
farthest_i = i;
}
}
counters[max_k]--;
counters[k]++;
labels[farthest_i] = k;
sample = data.ptr<float>(farthest_i);
for( j = 0; j < dims; j++ )
{
old_center[j] -= sample[j];
new_center[j] += sample[j];
}
}
for( k = 0; k < K; k++ )
{
float* center = centers.ptr<float>(k);
CV_Assert( counters[k] != 0 );
float scale = 1.f/counters[k];
for( j = 0; j < dims; j++ )
center[j] *= scale;
if( iter > 0 )
{
double dist = 0;
const float* old_center = old_centers.ptr<float>(k);
for( j = 0; j < dims; j++ )
{
double t = center[j] - old_center[j];
dist += t*t;
}
max_center_shift = std::max(max_center_shift, dist);
}
}
}
if( ++iter == MAX(criteria.maxCount, 2) || max_center_shift <= criteria.epsilon )
break;
// assign labels
Mat dists(1, N, CV_64F);
double* dist = dists.ptr<double>(0);
parallel_for_(Range(0, N),
KMeansDistanceComputer(dist, labels, data, centers));
compactness = 0;
for( i = 0; i < N; i++ )
{
compactness += dist[i];
}
}
if( compactness < best_compactness )
{
best_compactness = compactness;
if( _centers.needed() )
centers.copyTo(_centers);
_labels.copyTo(best_labels);
}
}
return best_compactness;
}

361
modules/core/src/lpsolver.cpp Normal file
View File

@@ -0,0 +1,361 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
#include <climits>
#include <algorithm>
#include <cstdarg>
#define dprintf(x)
#define print_matrix(x)
namespace cv{
using std::vector;
#ifdef ALEX_DEBUG
static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::vector<int> N,const std::vector<int> B){
printf("\tprint simplex state\n");
printf("v=%g\n",v);
printf("here c goes\n");
print_matrix(c);
printf("non-basic: ");
print(Mat(N));
printf("\n");
printf("here b goes\n");
print_matrix(b);
printf("basic: ");
print(Mat(B));
printf("\n");
}
#else
#define print_simplex_state(c,b,v,N,B)
#endif
/**Due to technical considerations, the format of input b and c is somewhat special:
*both b and c should be one column bigger than corresponding b and c of linear problem and the leftmost column will be used internally
by this procedure - it should not be cleaned before the call to procedure and may contain mess after
it also initializes N and B and does not make any assumptions about their init values
* @return SOLVELP_UNFEASIBLE if problem is unfeasible, 0 if feasible.
*/
static int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow);
static inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,int leaving_index,
int entering_index,vector<unsigned int>& indexToRow);
/**@return SOLVELP_UNBOUNDED means the problem is unbdd, SOLVELP_MULTI means multiple solutions, SOLVELP_SINGLE means one solution.
*/
static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow);
static void swap_columns(Mat_<double>& A,int col1,int col2);
#define SWAP(type,a,b) {type tmp=(a);(a)=(b);(b)=tmp;}
//return codes:-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm)
int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
dprintf(("call to solveLP\n"));
//sanity check (size, type, no. of channels)
CV_Assert(Func.type()==CV_64FC1 || Func.type()==CV_32FC1);
CV_Assert(Constr.type()==CV_64FC1 || Constr.type()==CV_32FC1);
CV_Assert((Func.rows==1 && (Constr.cols-Func.cols==1))||
(Func.cols==1 && (Constr.cols-Func.rows==1)));
//copy arguments for we will shall modify them
Mat_<double> bigC=Mat_<double>(1,(Func.rows==1?Func.cols:Func.rows)+1),
bigB=Mat_<double>(Constr.rows,Constr.cols+1);
if(Func.rows==1){
Func.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1);
}else{
Mat FuncT=Func.t();
FuncT.convertTo(bigC.colRange(1,bigC.cols),CV_64FC1);
}
Constr.convertTo(bigB.colRange(1,bigB.cols),CV_64FC1);
double v=0;
vector<int> N,B;
vector<unsigned int> indexToRow;
if(initialize_simplex(bigC,bigB,v,N,B,indexToRow)==SOLVELP_UNFEASIBLE){
return SOLVELP_UNFEASIBLE;
}
Mat_<double> c=bigC.colRange(1,bigC.cols),
b=bigB.colRange(1,bigB.cols);
int res=0;
if((res=inner_simplex(c,b,v,N,B,indexToRow))==SOLVELP_UNBOUNDED){
return SOLVELP_UNBOUNDED;
}
//return the optimal solution
z.create(c.cols,1,CV_64FC1);
MatIterator_<double> it=z.begin<double>();
unsigned int nsize = (unsigned int)N.size();
for(int i=1;i<=c.cols;i++,it++){
if(indexToRow[i]<nsize){
*it=0;
}else{
*it=b.at<double>(indexToRow[i]-nsize,b.cols-1);
}
}
return res;
}
static int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow){
N.resize(c.cols);
N[0]=0;
for (std::vector<int>::iterator it = N.begin()+1 ; it != N.end(); ++it){
*it=it[-1]+1;
}
B.resize(b.rows);
B[0]=(int)N.size();
for (std::vector<int>::iterator it = B.begin()+1 ; it != B.end(); ++it){
*it=it[-1]+1;
}
indexToRow.resize(c.cols+b.rows);
indexToRow[0]=0;
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
*it=it[-1]+1;
}
v=0;
int k=0;
{
double min=DBL_MAX;
for(int i=0;i<b.rows;i++){
if(b(i,b.cols-1)<min){
min=b(i,b.cols-1);
k=i;
}
}
}
if(b(k,b.cols-1)>=0){
N.erase(N.begin());
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
--(*it);
}
return 0;
}
Mat_<double> old_c=c.clone();
c=0;
c(0,0)=-1;
for(int i=0;i<b.rows;i++){
b(i,0)=-1;
}
print_simplex_state(c,b,v,N,B);
dprintf(("\tWE MAKE PIVOT\n"));
pivot(c,b,v,N,B,k,0,indexToRow);
print_simplex_state(c,b,v,N,B);
inner_simplex(c,b,v,N,B,indexToRow);
dprintf(("\tAFTER INNER_SIMPLEX\n"));
print_simplex_state(c,b,v,N,B);
unsigned int nsize = (unsigned int)N.size();
if(indexToRow[0]>=nsize){
int iterator_offset=indexToRow[0]-nsize;
if(b(iterator_offset,b.cols-1)>0){
return SOLVELP_UNFEASIBLE;
}
pivot(c,b,v,N,B,iterator_offset,0,indexToRow);
}
vector<int>::iterator iterator;
{
int iterator_offset=indexToRow[0];
iterator=N.begin()+iterator_offset;
std::iter_swap(iterator,N.begin());
SWAP(int,indexToRow[*iterator],indexToRow[0]);
swap_columns(c,iterator_offset,0);
swap_columns(b,iterator_offset,0);
}
dprintf(("after swaps\n"));
print_simplex_state(c,b,v,N,B);
//start from 1, because we ignore x_0
c=0;
v=0;
for(int I=1;I<old_c.cols;I++){
if(indexToRow[I]<nsize){
dprintf(("I=%d from nonbasic\n",I));
int iterator_offset=indexToRow[I];
c(0,iterator_offset)+=old_c(0,I);
print_matrix(c);
}else{
dprintf(("I=%d from basic\n",I));
int iterator_offset=indexToRow[I]-nsize;
c-=old_c(0,I)*b.row(iterator_offset).colRange(0,b.cols-1);
v+=old_c(0,I)*b(iterator_offset,b.cols-1);
print_matrix(c);
}
}
dprintf(("after restore\n"));
print_simplex_state(c,b,v,N,B);
N.erase(N.begin());
for (std::vector<unsigned int>::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){
--(*it);
}
return 0;
}
static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow){
int count=0;
for(;;){
dprintf(("iteration #%d\n",count));
count++;
static MatIterator_<double> pos_ptr;
int e=-1,pos_ctr=0,min_var=INT_MAX;
bool all_nonzero=true;
for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){
if(*pos_ptr==0){
all_nonzero=false;
}
if(*pos_ptr>0){
if(N[pos_ctr]<min_var){
e=pos_ctr;
min_var=N[pos_ctr];
}
}
}
if(e==-1){
dprintf(("hello from e==-1\n"));
print_matrix(c);
if(all_nonzero==true){
return SOLVELP_SINGLE;
}else{
return SOLVELP_MULTI;
}
}
int l=-1;
min_var=INT_MAX;
double min=DBL_MAX;
int row_it=0;
MatIterator_<double> min_row_ptr=b.begin();
for(MatIterator_<double> it=b.begin();it!=b.end();it+=b.cols,row_it++){
double myite=0;
//check constraints, select the tightest one, reinforcing Bland's rule
if((myite=it[e])>0){
double val=it[b.cols-1]/myite;
if(val<min || (val==min && B[row_it]<min_var)){
min_var=B[row_it];
min_row_ptr=it;
min=val;
l=row_it;
}
}
}
if(l==-1){
return SOLVELP_UNBOUNDED;
}
dprintf(("the tightest constraint is in row %d with %g\n",l,min));
pivot(c,b,v,N,B,l,e,indexToRow);
dprintf(("objective, v=%g\n",v));
print_matrix(c);
dprintf(("constraints\n"));
print_matrix(b);
dprintf(("non-basic: "));
print_matrix(Mat(N));
dprintf(("basic: "));
print_matrix(Mat(B));
}
}
static inline void pivot(Mat_<double>& c,Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,
int leaving_index,int entering_index,vector<unsigned int>& indexToRow){
double Coef=b(leaving_index,entering_index);
for(int i=0;i<b.cols;i++){
if(i==entering_index){
b(leaving_index,i)=1/Coef;
}else{
b(leaving_index,i)/=Coef;
}
}
for(int i=0;i<b.rows;i++){
if(i!=leaving_index){
double coef=b(i,entering_index);
for(int j=0;j<b.cols;j++){
if(j==entering_index){
b(i,j)=-coef*b(leaving_index,j);
}else{
b(i,j)-=(coef*b(leaving_index,j));
}
}
}
}
//objective function
Coef=c(0,entering_index);
for(int i=0;i<(b.cols-1);i++){
if(i==entering_index){
c(0,i)=-Coef*b(leaving_index,i);
}else{
c(0,i)-=Coef*b(leaving_index,i);
}
}
dprintf(("v was %g\n",v));
v+=Coef*b(leaving_index,b.cols-1);
SWAP(int,N[entering_index],B[leaving_index]);
SWAP(int,indexToRow[N[entering_index]],indexToRow[B[leaving_index]]);
}
static inline void swap_columns(Mat_<double>& A,int col1,int col2){
for(int i=0;i<A.rows;i++){
double tmp=A(i,col1);
A(i,col1)=A(i,col2);
A(i,col2)=tmp;
}
}
}

335
modules/core/src/matmul.cpp
View File

@@ -2959,341 +2959,6 @@ double Mat::dot(InputArray _mat) const
return r;
}
/****************************************************************************************\
* PCA *
\****************************************************************************************/
PCA::PCA() {}
PCA::PCA(InputArray data, InputArray _mean, int flags, int maxComponents)
{
operator()(data, _mean, flags, maxComponents);
}
PCA::PCA(InputArray data, InputArray _mean, int flags, double retainedVariance)
{
operator()(data, _mean, flags, retainedVariance);
}
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents)
{
Mat data = _data.getMat(), _mean = __mean.getMat();
int covar_flags = CV_COVAR_SCALE;
int i, len, in_count;
Size mean_sz;
CV_Assert( data.channels() == 1 );
if( flags & CV_PCA_DATA_AS_COL )
{
len = data.rows;
in_count = data.cols;
covar_flags |= CV_COVAR_COLS;
mean_sz = Size(1, len);
}
else
{
len = data.cols;
in_count = data.rows;
covar_flags |= CV_COVAR_ROWS;
mean_sz = Size(len, 1);
}
int count = std::min(len, in_count), out_count = count;
if( maxComponents > 0 )
out_count = std::min(count, maxComponents);
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
if( len <= in_count )
covar_flags |= CV_COVAR_NORMAL;
int ctype = std::max(CV_32F, data.depth());
mean.create( mean_sz, ctype );
Mat covar( count, count, ctype );
if( !_mean.empty() )
{
CV_Assert( _mean.size() == mean_sz );
_mean.convertTo(mean, ctype);
covar_flags |= CV_COVAR_USE_AVG;
}
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
eigen( covar, eigenvalues, eigenvectors );
if( !(covar_flags & CV_COVAR_NORMAL) )
{
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
Mat evects1(count, len, ctype);
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
eigenvectors = evects1;
// normalize eigenvectors
for( i = 0; i < out_count; i++ )
{
Mat vec = eigenvectors.row(i);
normalize(vec, vec);
}
}
if( count > out_count )
{
// use clone() to physically copy the data and thus deallocate the original matrices
eigenvalues = eigenvalues.rowRange(0,out_count).clone();
eigenvectors = eigenvectors.rowRange(0,out_count).clone();
}
return *this;
}
void PCA::write(FileStorage& fs ) const
{
CV_Assert( fs.isOpened() );
fs << "name" << "PCA";
fs << "vectors" << eigenvectors;
fs << "values" << eigenvalues;
fs << "mean" << mean;
}
void PCA::read(const FileNode& fs)
{
CV_Assert( !fs.empty() );
String name = (String)fs["name"];
CV_Assert( name == "PCA" );
cv::read(fs["vectors"], eigenvectors);
cv::read(fs["values"], eigenvalues);
cv::read(fs["mean"], mean);
}
template <typename T>
int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance)
{
CV_DbgAssert( eigenvalues.type() == DataType<T>::type );
Mat g(eigenvalues.size(), DataType<T>::type);
for(int ig = 0; ig < g.rows; ig++)
{
g.at<T>(ig, 0) = 0;
for(int im = 0; im <= ig; im++)
{
g.at<T>(ig,0) += eigenvalues.at<T>(im,0);
}
}
int L;
for(L = 0; L < eigenvalues.rows; L++)
{
double energy = g.at<T>(L, 0) / g.at<T>(g.rows - 1, 0);
if(energy > retainedVariance)
break;
}
L = std::max(2, L);
return L;
}
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, double retainedVariance)
{
Mat data = _data.getMat(), _mean = __mean.getMat();
int covar_flags = CV_COVAR_SCALE;
int i, len, in_count;
Size mean_sz;
CV_Assert( data.channels() == 1 );
if( flags & CV_PCA_DATA_AS_COL )
{
len = data.rows;
in_count = data.cols;
covar_flags |= CV_COVAR_COLS;
mean_sz = Size(1, len);
}
else
{
len = data.cols;
in_count = data.rows;
covar_flags |= CV_COVAR_ROWS;
mean_sz = Size(len, 1);
}
CV_Assert( retainedVariance > 0 && retainedVariance <= 1 );
int count = std::min(len, in_count);
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
if( len <= in_count )
covar_flags |= CV_COVAR_NORMAL;
int ctype = std::max(CV_32F, data.depth());
mean.create( mean_sz, ctype );
Mat covar( count, count, ctype );
if( !_mean.empty() )
{
CV_Assert( _mean.size() == mean_sz );
_mean.convertTo(mean, ctype);
}
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
eigen( covar, eigenvalues, eigenvectors );
if( !(covar_flags & CV_COVAR_NORMAL) )
{
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
Mat evects1(count, len, ctype);
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
eigenvectors = evects1;
// normalize all eigenvectors
for( i = 0; i < eigenvectors.rows; i++ )
{
Mat vec = eigenvectors.row(i);
normalize(vec, vec);
}
}
// compute the cumulative energy content for each eigenvector
int L;
if (ctype == CV_32F)
L = computeCumulativeEnergy<float>(eigenvalues, retainedVariance);
else
L = computeCumulativeEnergy<double>(eigenvalues, retainedVariance);
// use clone() to physically copy the data and thus deallocate the original matrices
eigenvalues = eigenvalues.rowRange(0,L).clone();
eigenvectors = eigenvectors.rowRange(0,L).clone();
return *this;
}
void PCA::project(InputArray _data, OutputArray result) const
{
Mat data = _data.getMat();
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
((mean.rows == 1 && mean.cols == data.cols) || (mean.cols == 1 && mean.rows == data.rows)));
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
int ctype = mean.type();
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
if( mean.rows == 1 )
gemm( tmp_data, eigenvectors, 1, Mat(), 0, result, GEMM_2_T );
else
gemm( eigenvectors, tmp_data, 1, Mat(), 0, result, 0 );
}
Mat PCA::project(InputArray data) const
{
Mat result;
project(data, result);
return result;
}
void PCA::backProject(InputArray _data, OutputArray result) const
{
Mat data = _data.getMat();
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
((mean.rows == 1 && eigenvectors.rows == data.cols) ||
(mean.cols == 1 && eigenvectors.rows == data.rows)));
Mat tmp_data, tmp_mean;
data.convertTo(tmp_data, mean.type());
if( mean.rows == 1 )
{
tmp_mean = repeat(mean, data.rows, 1);
gemm( tmp_data, eigenvectors, 1, tmp_mean, 1, result, 0 );
}
else
{
tmp_mean = repeat(mean, 1, data.cols);
gemm( eigenvectors, tmp_data, 1, tmp_mean, 1, result, GEMM_1_T );
}
}
Mat PCA::backProject(InputArray data) const
{
Mat result;
backProject(data, result);
return result;
}
}
void cv::PCACompute(InputArray data, InputOutputArray mean,
OutputArray eigenvectors, int maxComponents)
{
PCA pca;
pca(data, mean, 0, maxComponents);
pca.mean.copyTo(mean);
pca.eigenvectors.copyTo(eigenvectors);
}
void cv::PCACompute(InputArray data, InputOutputArray mean,
OutputArray eigenvectors, double retainedVariance)
{
PCA pca;
pca(data, mean, 0, retainedVariance);
pca.mean.copyTo(mean);
pca.eigenvectors.copyTo(eigenvectors);
}
void cv::PCAProject(InputArray data, InputArray mean,
InputArray eigenvectors, OutputArray result)
{
PCA pca;
pca.mean = mean.getMat();
pca.eigenvectors = eigenvectors.getMat();
pca.project(data, result);
}
void cv::PCABackProject(InputArray data, InputArray mean,
InputArray eigenvectors, OutputArray result)
{
PCA pca;
pca.mean = mean.getMat();
pca.eigenvectors = eigenvectors.getMat();
pca.backProject(data, result);
}
/****************************************************************************************\

414
modules/core/src/matrix.cpp
View File

@@ -3971,420 +3971,6 @@ void cv::sortIdx( InputArray _src, OutputArray _dst, int flags )
}
////////////////////////////////////////// kmeans ////////////////////////////////////////////
namespace cv
{
static void generateRandomCenter(const std::vector<Vec2f>& box, float* center, RNG& rng)
{
size_t j, dims = box.size();
float margin = 1.f/dims;
for( j = 0; j < dims; j++ )
center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
}
class KMeansPPDistanceComputer : public ParallelLoopBody
{
public:
KMeansPPDistanceComputer( float *_tdist2,
const float *_data,
const float *_dist,
int _dims,
size_t _step,
size_t _stepci )
: tdist2(_tdist2),
data(_data),
dist(_dist),
dims(_dims),
step(_step),
stepci(_stepci) { }
void operator()( const cv::Range& range ) const
{
const int begin = range.start;
const int end = range.end;
for ( int i = begin; i<end; i++ )
{
tdist2[i] = std::min(normL2Sqr_(data + step*i, data + stepci, dims), dist[i]);
}
}
private:
KMeansPPDistanceComputer& operator=(const KMeansPPDistanceComputer&); // to quiet MSVC
float *tdist2;
const float *data;
const float *dist;
const int dims;
const size_t step;
const size_t stepci;
};
/*
k-means center initialization using the following algorithm:
Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding
*/
static void generateCentersPP(const Mat& _data, Mat& _out_centers,
int K, RNG& rng, int trials)
{
int i, j, k, dims = _data.cols, N = _data.rows;
const float* data = _data.ptr<float>(0);
size_t step = _data.step/sizeof(data[0]);
std::vector<int> _centers(K);
int* centers = &_centers[0];
std::vector<float> _dist(N*3);
float* dist = &_dist[0], *tdist = dist + N, *tdist2 = tdist + N;
double sum0 = 0;
centers[0] = (unsigned)rng % N;
for( i = 0; i < N; i++ )
{
dist[i] = normL2Sqr_(data + step*i, data + step*centers[0], dims);
sum0 += dist[i];
}
for( k = 1; k < K; k++ )
{
double bestSum = DBL_MAX;
int bestCenter = -1;
for( j = 0; j < trials; j++ )
{
double p = (double)rng*sum0, s = 0;
for( i = 0; i < N-1; i++ )
if( (p -= dist[i]) <= 0 )
break;
int ci = i;
parallel_for_(Range(0, N),
KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
for( i = 0; i < N; i++ )
{
s += tdist2[i];
}
if( s < bestSum )
{
bestSum = s;
bestCenter = ci;
std::swap(tdist, tdist2);
}
}
centers[k] = bestCenter;
sum0 = bestSum;
std::swap(dist, tdist);
}
for( k = 0; k < K; k++ )
{
const float* src = data + step*centers[k];
float* dst = _out_centers.ptr<float>(k);
for( j = 0; j < dims; j++ )
dst[j] = src[j];
}
}
class KMeansDistanceComputer : public ParallelLoopBody
{
public:
KMeansDistanceComputer( double *_distances,
int *_labels,
const Mat& _data,
const Mat& _centers )
: distances(_distances),
labels(_labels),
data(_data),
centers(_centers)
{
}
void operator()( const Range& range ) const
{
const int begin = range.start;
const int end = range.end;
const int K = centers.rows;
const int dims = centers.cols;
const float *sample;
for( int i = begin; i<end; ++i)
{
sample = data.ptr<float>(i);
int k_best = 0;
double min_dist = DBL_MAX;
for( int k = 0; k < K; k++ )
{
const float* center = centers.ptr<float>(k);
const double dist = normL2Sqr_(sample, center, dims);
if( min_dist > dist )
{
min_dist = dist;
k_best = k;
}
}
distances[i] = min_dist;
labels[i] = k_best;
}
}
private:
KMeansDistanceComputer& operator=(const KMeansDistanceComputer&); // to quiet MSVC
double *distances;
int *labels;
const Mat& data;
const Mat& centers;
};
}
double cv::kmeans( InputArray _data, int K,
InputOutputArray _bestLabels,
TermCriteria criteria, int attempts,
int flags, OutputArray _centers )
{
const int SPP_TRIALS = 3;
Mat data0 = _data.getMat();
bool isrow = data0.rows == 1 && data0.channels() > 1;
int N = !isrow ? data0.rows : data0.cols;
int dims = (!isrow ? data0.cols : 1)*data0.channels();
int type = data0.depth();
attempts = std::max(attempts, 1);
CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
CV_Assert( N >= K );
Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));
_bestLabels.create(N, 1, CV_32S, -1, true);
Mat _labels, best_labels = _bestLabels.getMat();
if( flags & CV_KMEANS_USE_INITIAL_LABELS )
{
CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
best_labels.cols*best_labels.rows == N &&
best_labels.type() == CV_32S &&
best_labels.isContinuous());
best_labels.copyTo(_labels);
}
else
{
if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
best_labels.cols*best_labels.rows == N &&
best_labels.type() == CV_32S &&
best_labels.isContinuous()))
best_labels.create(N, 1, CV_32S);
_labels.create(best_labels.size(), best_labels.type());
}
int* labels = _labels.ptr<int>();
Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
std::vector<int> counters(K);
std::vector<Vec2f> _box(dims);
Vec2f* box = &_box[0];
double best_compactness = DBL_MAX, compactness = 0;
RNG& rng = theRNG();
int a, iter, i, j, k;
if( criteria.type & TermCriteria::EPS )
criteria.epsilon = std::max(criteria.epsilon, 0.);
else
criteria.epsilon = FLT_EPSILON;
criteria.epsilon *= criteria.epsilon;
if( criteria.type & TermCriteria::COUNT )
criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
else
criteria.maxCount = 100;
if( K == 1 )
{
attempts = 1;
criteria.maxCount = 2;
}
const float* sample = data.ptr<float>(0);
for( j = 0; j < dims; j++ )
box[j] = Vec2f(sample[j], sample[j]);
for( i = 1; i < N; i++ )
{
sample = data.ptr<float>(i);
for( j = 0; j < dims; j++ )
{
float v = sample[j];
box[j][0] = std::min(box[j][0], v);
box[j][1] = std::max(box[j][1], v);
}
}
for( a = 0; a < attempts; a++ )
{
double max_center_shift = DBL_MAX;
for( iter = 0;; )
{
swap(centers, old_centers);
if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
{
if( flags & KMEANS_PP_CENTERS )
generateCentersPP(data, centers, K, rng, SPP_TRIALS);
else
{
for( k = 0; k < K; k++ )
generateRandomCenter(_box, centers.ptr<float>(k), rng);
}
}
else
{
if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
{
for( i = 0; i < N; i++ )
CV_Assert( (unsigned)labels[i] < (unsigned)K );
}
// compute centers
centers = Scalar(0);
for( k = 0; k < K; k++ )
counters[k] = 0;
for( i = 0; i < N; i++ )
{
sample = data.ptr<float>(i);
k = labels[i];
float* center = centers.ptr<float>(k);
j=0;
#if CV_ENABLE_UNROLLED
for(; j <= dims - 4; j += 4 )
{
float t0 = center[j] + sample[j];
float t1 = center[j+1] + sample[j+1];
center[j] = t0;
center[j+1] = t1;
t0 = center[j+2] + sample[j+2];
t1 = center[j+3] + sample[j+3];
center[j+2] = t0;
center[j+3] = t1;
}
#endif
for( ; j < dims; j++ )
center[j] += sample[j];
counters[k]++;
}
if( iter > 0 )
max_center_shift = 0;
for( k = 0; k < K; k++ )
{
if( counters[k] != 0 )
continue;
// if some cluster appeared to be empty then:
// 1. find the biggest cluster
// 2. find the farthest from the center point in the biggest cluster
// 3. exclude the farthest point from the biggest cluster and form a new 1-point cluster.
int max_k = 0;
for( int k1 = 1; k1 < K; k1++ )
{
if( counters[max_k] < counters[k1] )
max_k = k1;
}
double max_dist = 0;
int farthest_i = -1;
float* new_center = centers.ptr<float>(k);
float* old_center = centers.ptr<float>(max_k);
float* _old_center = temp.ptr<float>(); // normalized
float scale = 1.f/counters[max_k];
for( j = 0; j < dims; j++ )
_old_center[j] = old_center[j]*scale;
for( i = 0; i < N; i++ )
{
if( labels[i] != max_k )
continue;
sample = data.ptr<float>(i);
double dist = normL2Sqr_(sample, _old_center, dims);
if( max_dist <= dist )
{
max_dist = dist;
farthest_i = i;
}
}
counters[max_k]--;
counters[k]++;
labels[farthest_i] = k;
sample = data.ptr<float>(farthest_i);
for( j = 0; j < dims; j++ )
{
old_center[j] -= sample[j];
new_center[j] += sample[j];
}
}
for( k = 0; k < K; k++ )
{
float* center = centers.ptr<float>(k);
CV_Assert( counters[k] != 0 );
float scale = 1.f/counters[k];
for( j = 0; j < dims; j++ )
center[j] *= scale;
if( iter > 0 )
{
double dist = 0;
const float* old_center = old_centers.ptr<float>(k);
for( j = 0; j < dims; j++ )
{
double t = center[j] - old_center[j];
dist += t*t;
}
max_center_shift = std::max(max_center_shift, dist);
}
}
}
if( ++iter == MAX(criteria.maxCount, 2) || max_center_shift <= criteria.epsilon )
break;
// assign labels
Mat dists(1, N, CV_64F);
double* dist = dists.ptr<double>(0);
parallel_for_(Range(0, N),
KMeansDistanceComputer(dist, labels, data, centers));
compactness = 0;
for( i = 0; i < N; i++ )
{
compactness += dist[i];
}
}
if( compactness < best_compactness )
{
best_compactness = compactness;
if( _centers.needed() )
centers.copyTo(_centers);
_labels.copyTo(best_labels);
}
}
return best_compactness;
}
CV_IMPL void cvSetIdentity( CvArr* arr, CvScalar value )
{
cv::Mat m = cv::cvarrToMat(arr);

384
modules/core/src/pca.cpp Normal file
View File

@@ -0,0 +1,384 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "precomp.hpp"
/****************************************************************************************\
* PCA *
\****************************************************************************************/
namespace cv
{
PCA::PCA() {}
PCA::PCA(InputArray data, InputArray _mean, int flags, int maxComponents)
{
operator()(data, _mean, flags, maxComponents);
}
PCA::PCA(InputArray data, InputArray _mean, int flags, double retainedVariance)
{
operator()(data, _mean, flags, retainedVariance);
}
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents)
{
Mat data = _data.getMat(), _mean = __mean.getMat();
int covar_flags = CV_COVAR_SCALE;
int i, len, in_count;
Size mean_sz;
CV_Assert( data.channels() == 1 );
if( flags & CV_PCA_DATA_AS_COL )
{
len = data.rows;
in_count = data.cols;
covar_flags |= CV_COVAR_COLS;
mean_sz = Size(1, len);
}
else
{
len = data.cols;
in_count = data.rows;
covar_flags |= CV_COVAR_ROWS;
mean_sz = Size(len, 1);
}
int count = std::min(len, in_count), out_count = count;
if( maxComponents > 0 )
out_count = std::min(count, maxComponents);
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
if( len <= in_count )
covar_flags |= CV_COVAR_NORMAL;
int ctype = std::max(CV_32F, data.depth());
mean.create( mean_sz, ctype );
Mat covar( count, count, ctype );
if( !_mean.empty() )
{
CV_Assert( _mean.size() == mean_sz );
_mean.convertTo(mean, ctype);
covar_flags |= CV_COVAR_USE_AVG;
}
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
eigen( covar, eigenvalues, eigenvectors );
if( !(covar_flags & CV_COVAR_NORMAL) )
{
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
Mat evects1(count, len, ctype);
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
eigenvectors = evects1;
// normalize eigenvectors
for( i = 0; i < out_count; i++ )
{
Mat vec = eigenvectors.row(i);
normalize(vec, vec);
}
}
if( count > out_count )
{
// use clone() to physically copy the data and thus deallocate the original matrices
eigenvalues = eigenvalues.rowRange(0,out_count).clone();
eigenvectors = eigenvectors.rowRange(0,out_count).clone();
}
return *this;
}
void PCA::write(FileStorage& fs ) const
{
CV_Assert( fs.isOpened() );
fs << "name" << "PCA";
fs << "vectors" << eigenvectors;
fs << "values" << eigenvalues;
fs << "mean" << mean;
}
void PCA::read(const FileNode& fs)
{
CV_Assert( !fs.empty() );
String name = (String)fs["name"];
CV_Assert( name == "PCA" );
cv::read(fs["vectors"], eigenvectors);
cv::read(fs["values"], eigenvalues);
cv::read(fs["mean"], mean);
}
template <typename T>
int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance)
{
CV_DbgAssert( eigenvalues.type() == DataType<T>::type );
Mat g(eigenvalues.size(), DataType<T>::type);
for(int ig = 0; ig < g.rows; ig++)
{
g.at<T>(ig, 0) = 0;
for(int im = 0; im <= ig; im++)
{
g.at<T>(ig,0) += eigenvalues.at<T>(im,0);
}
}
int L;
for(L = 0; L < eigenvalues.rows; L++)
{
double energy = g.at<T>(L, 0) / g.at<T>(g.rows - 1, 0);
if(energy > retainedVariance)
break;
}
L = std::max(2, L);
return L;
}
PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, double retainedVariance)
{
Mat data = _data.getMat(), _mean = __mean.getMat();
int covar_flags = CV_COVAR_SCALE;
int i, len, in_count;
Size mean_sz;
CV_Assert( data.channels() == 1 );
if( flags & CV_PCA_DATA_AS_COL )
{
len = data.rows;
in_count = data.cols;
covar_flags |= CV_COVAR_COLS;
mean_sz = Size(1, len);
}
else
{
len = data.cols;
in_count = data.rows;
covar_flags |= CV_COVAR_ROWS;
mean_sz = Size(len, 1);
}
CV_Assert( retainedVariance > 0 && retainedVariance <= 1 );
int count = std::min(len, in_count);
// "scrambled" way to compute PCA (when cols(A)>rows(A)):
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
if( len <= in_count )
covar_flags |= CV_COVAR_NORMAL;
int ctype = std::max(CV_32F, data.depth());
mean.create( mean_sz, ctype );
Mat covar( count, count, ctype );
if( !_mean.empty() )
{
CV_Assert( _mean.size() == mean_sz );
_mean.convertTo(mean, ctype);
}
calcCovarMatrix( data, covar, mean, covar_flags, ctype );
eigen( covar, eigenvalues, eigenvectors );
if( !(covar_flags & CV_COVAR_NORMAL) )
{
// CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
// CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
Mat evects1(count, len, ctype);
gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
(flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
eigenvectors = evects1;
// normalize all eigenvectors
for( i = 0; i < eigenvectors.rows; i++ )
{
Mat vec = eigenvectors.row(i);
normalize(vec, vec);
}
}
// compute the cumulative energy content for each eigenvector
int L;
if (ctype == CV_32F)
L = computeCumulativeEnergy<float>(eigenvalues, retainedVariance);
else
L = computeCumulativeEnergy<double>(eigenvalues, retainedVariance);
// use clone() to physically copy the data and thus deallocate the original matrices
eigenvalues = eigenvalues.rowRange(0,L).clone();
eigenvectors = eigenvectors.rowRange(0,L).clone();
return *this;
}
void PCA::project(InputArray _data, OutputArray result) const
{
Mat data = _data.getMat();
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
((mean.rows == 1 && mean.cols == data.cols) || (mean.cols == 1 && mean.rows == data.rows)));
Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
int ctype = mean.type();
if( data.type() != ctype || tmp_mean.data == mean.data )
{
data.convertTo( tmp_data, ctype );
subtract( tmp_data, tmp_mean, tmp_data );
}
else
{
subtract( data, tmp_mean, tmp_mean );
tmp_data = tmp_mean;
}
if( mean.rows == 1 )
gemm( tmp_data, eigenvectors, 1, Mat(), 0, result, GEMM_2_T );
else
gemm( eigenvectors, tmp_data, 1, Mat(), 0, result, 0 );
}
Mat PCA::project(InputArray data) const
{
Mat result;
project(data, result);
return result;
}
void PCA::backProject(InputArray _data, OutputArray result) const
{
Mat data = _data.getMat();
CV_Assert( !mean.empty() && !eigenvectors.empty() &&
((mean.rows == 1 && eigenvectors.rows == data.cols) ||
(mean.cols == 1 && eigenvectors.rows == data.rows)));
Mat tmp_data, tmp_mean;
data.convertTo(tmp_data, mean.type());
if( mean.rows == 1 )
{
tmp_mean = repeat(mean, data.rows, 1);
gemm( tmp_data, eigenvectors, 1, tmp_mean, 1, result, 0 );
}
else
{
tmp_mean = repeat(mean, 1, data.cols);
gemm( eigenvectors, tmp_data, 1, tmp_mean, 1, result, GEMM_1_T );
}
}
Mat PCA::backProject(InputArray data) const
{
Mat result;
backProject(data, result);
return result;
}
}
void cv::PCACompute(InputArray data, InputOutputArray mean,
OutputArray eigenvectors, int maxComponents)
{
PCA pca;
pca(data, mean, 0, maxComponents);
pca.mean.copyTo(mean);
pca.eigenvectors.copyTo(eigenvectors);
}
void cv::PCACompute(InputArray data, InputOutputArray mean,
OutputArray eigenvectors, double retainedVariance)
{
PCA pca;
pca(data, mean, 0, retainedVariance);
pca.mean.copyTo(mean);
pca.eigenvectors.copyTo(eigenvectors);
}
void cv::PCAProject(InputArray data, InputArray mean,
InputArray eigenvectors, OutputArray result)
{
PCA pca;
pca.mean = mean.getMat();
pca.eigenvectors = eigenvectors.getMat();
pca.project(data, result);
}
void cv::PCABackProject(InputArray data, InputArray mean,
InputArray eigenvectors, OutputArray result)
{
PCA pca;
pca.mean = mean.getMat();
pca.eigenvectors = eigenvectors.getMat();
pca.backProject(data, result);
}

102
modules/core/test/test_conjugate_gradient.cpp Normal file
View File

@@ -0,0 +1,102 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "test_precomp.hpp"
#include <cstdlib>
static void mytest(cv::Ptr<cv::ConjGradSolver> solver,cv::Ptr<cv::MinProblemSolver::Function> ptr_F,cv::Mat& x,
cv::Mat& etalon_x,double etalon_res){
solver->setFunction(ptr_F);
//int ndim=MAX(step.cols,step.rows);
double res=solver->minimize(x);
std::cout<<"res:\n\t"<<res<<std::endl;
std::cout<<"x:\n\t"<<x<<std::endl;
std::cout<<"etalon_res:\n\t"<<etalon_res<<std::endl;
std::cout<<"etalon_x:\n\t"<<etalon_x<<std::endl;
double tol=solver->getTermCriteria().epsilon;
ASSERT_TRUE(std::abs(res-etalon_res)<tol);
/*for(cv::Mat_<double>::iterator it1=x.begin<double>(),it2=etalon_x.begin<double>();it1!=x.end<double>();it1++,it2++){
ASSERT_TRUE(std::abs((*it1)-(*it2))<tol);
}*/
std::cout<<"--------------------------\n";
}
class SphereF:public cv::MinProblemSolver::Function{
public:
double calc(const double* x)const{
return x[0]*x[0]+x[1]*x[1]+x[2]*x[2]+x[3]*x[3];
}
void getGradient(const double* x,double* grad){
for(int i=0;i<4;i++){
grad[i]=2*x[i];
}
}
};
class RosenbrockF:public cv::MinProblemSolver::Function{
double calc(const double* x)const{
return 100*(x[1]-x[0]*x[0])*(x[1]-x[0]*x[0])+(1-x[0])*(1-x[0]);
}
void getGradient(const double* x,double* grad){
grad[0]=-2*(1-x[0])-400*(x[1]-x[0]*x[0])*x[0];
grad[1]=200*(x[1]-x[0]*x[0]);
}
};
TEST(Optim_ConjGrad, regression_basic){
cv::Ptr<cv::ConjGradSolver> solver=cv::ConjGradSolver::create();
#if 1
{
cv::Ptr<cv::MinProblemSolver::Function> ptr_F(new SphereF());
cv::Mat x=(cv::Mat_<double>(4,1)<<50.0,10.0,1.0,-10.0),
etalon_x=(cv::Mat_<double>(1,4)<<0.0,0.0,0.0,0.0);
double etalon_res=0.0;
mytest(solver,ptr_F,x,etalon_x,etalon_res);
}
#endif
#if 1
{
cv::Ptr<cv::MinProblemSolver::Function> ptr_F(new RosenbrockF());
cv::Mat x=(cv::Mat_<double>(2,1)<<0.0,0.0),
etalon_x=(cv::Mat_<double>(2,1)<<1.0,1.0);
double etalon_res=0.0;
mytest(solver,ptr_F,x,etalon_x,etalon_res);
}
#endif
}

103
modules/core/test/test_downhill_simplex.cpp Normal file
View File

@@ -0,0 +1,103 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "test_precomp.hpp"
#include <cstdlib>
#include <cmath>
#include <algorithm>
static void mytest(cv::Ptr<cv::DownhillSolver> solver,cv::Ptr<cv::MinProblemSolver::Function> ptr_F,cv::Mat& x,cv::Mat& step,
cv::Mat& etalon_x,double etalon_res){
solver->setFunction(ptr_F);
int ndim=MAX(step.cols,step.rows);
solver->setInitStep(step);
cv::Mat settedStep;
solver->getInitStep(settedStep);
ASSERT_TRUE(settedStep.rows==1 && settedStep.cols==ndim);
ASSERT_TRUE(std::equal(step.begin<double>(),step.end<double>(),settedStep.begin<double>()));
std::cout<<"step setted:\n\t"<<step<<std::endl;
double res=solver->minimize(x);
std::cout<<"res:\n\t"<<res<<std::endl;
std::cout<<"x:\n\t"<<x<<std::endl;
std::cout<<"etalon_res:\n\t"<<etalon_res<<std::endl;
std::cout<<"etalon_x:\n\t"<<etalon_x<<std::endl;
double tol=solver->getTermCriteria().epsilon;
ASSERT_TRUE(std::abs(res-etalon_res)<tol);
/*for(cv::Mat_<double>::iterator it1=x.begin<double>(),it2=etalon_x.begin<double>();it1!=x.end<double>();it1++,it2++){
ASSERT_TRUE(std::abs((*it1)-(*it2))<tol);
}*/
std::cout<<"--------------------------\n";
}
class SphereF:public cv::MinProblemSolver::Function{
public:
double calc(const double* x)const{
return x[0]*x[0]+x[1]*x[1];
}
};
class RosenbrockF:public cv::MinProblemSolver::Function{
double calc(const double* x)const{
return 100*(x[1]-x[0]*x[0])*(x[1]-x[0]*x[0])+(1-x[0])*(1-x[0]);
}
};
TEST(Optim_Downhill, regression_basic){
cv::Ptr<cv::DownhillSolver> solver=cv::DownhillSolver::create();
#if 1
{
cv::Ptr<cv::MinProblemSolver::Function> ptr_F = cv::makePtr<SphereF>();
cv::Mat x=(cv::Mat_<double>(1,2)<<1.0,1.0),
step=(cv::Mat_<double>(2,1)<<-0.5,-0.5),
etalon_x=(cv::Mat_<double>(1,2)<<-0.0,0.0);
double etalon_res=0.0;
mytest(solver,ptr_F,x,step,etalon_x,etalon_res);
}
#endif
#if 1
{
cv::Ptr<cv::MinProblemSolver::Function> ptr_F = cv::makePtr<RosenbrockF>();
cv::Mat x=(cv::Mat_<double>(2,1)<<0.0,0.0),
step=(cv::Mat_<double>(2,1)<<0.5,+0.5),
etalon_x=(cv::Mat_<double>(2,1)<<1.0,1.0);
double etalon_res=0.0;
mytest(solver,ptr_F,x,step,etalon_x,etalon_res);
}
#endif
}

141
modules/core/test/test_lpsolver.cpp Normal file
View File

@@ -0,0 +1,141 @@
/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of the copyright holders may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the OpenCV Foundation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include "test_precomp.hpp"
#include <iostream>
TEST(Optim_LpSolver, regression_basic){
cv::Mat A,B,z,etalon_z;
#if 1
//cormen's example #1
A=(cv::Mat_<double>(3,1)<<3,1,2);
B=(cv::Mat_<double>(3,4)<<1,1,3,30,2,2,5,24,4,1,2,36);
std::cout<<"here A goes\n"<<A<<"\n";
cv::solveLP(A,B,z);
std::cout<<"here z goes\n"<<z<<"\n";
etalon_z=(cv::Mat_<double>(3,1)<<8,4,0);
ASSERT_EQ(cv::countNonZero(z!=etalon_z),0);
#endif
#if 1
//cormen's example #2
A=(cv::Mat_<double>(1,2)<<18,12.5);
B=(cv::Mat_<double>(3,3)<<1,1,20,1,0,20,0,1,16);
std::cout<<"here A goes\n"<<A<<"\n";
cv::solveLP(A,B,z);
std::cout<<"here z goes\n"<<z<<"\n";
etalon_z=(cv::Mat_<double>(2,1)<<20,0);
ASSERT_EQ(cv::countNonZero(z!=etalon_z),0);
#endif
#if 1
//cormen's example #3
A=(cv::Mat_<double>(1,2)<<5,-3);
B=(cv::Mat_<double>(2,3)<<1,-1,1,2,1,2);
std::cout<<"here A goes\n"<<A<<"\n";
cv::solveLP(A,B,z);
std::cout<<"here z goes\n"<<z<<"\n";
etalon_z=(cv::Mat_<double>(2,1)<<1,0);
ASSERT_EQ(cv::countNonZero(z!=etalon_z),0);
#endif
}
TEST(Optim_LpSolver, regression_init_unfeasible){
cv::Mat A,B,z,etalon_z;
#if 1
//cormen's example #4 - unfeasible
A=(cv::Mat_<double>(1,3)<<-1,-1,-1);
B=(cv::Mat_<double>(2,4)<<-2,-7.5,-3,-10000,-20,-5,-10,-30000);
std::cout<<"here A goes\n"<<A<<"\n";
cv::solveLP(A,B,z);
std::cout<<"here z goes\n"<<z<<"\n";
etalon_z=(cv::Mat_<double>(3,1)<<1250,1000,0);
ASSERT_EQ(cv::countNonZero(z!=etalon_z),0);
#endif
}
TEST(Optim_LpSolver, regression_absolutely_unfeasible){
cv::Mat A,B,z,etalon_z;
#if 1
//trivial absolutely unfeasible example
A=(cv::Mat_<double>(1,1)<<1);
B=(cv::Mat_<double>(2,2)<<1,-1);
std::cout<<"here A goes\n"<<A<<"\n";
int res=cv::solveLP(A,B,z);
ASSERT_EQ(res,-1);
#endif
}
TEST(Optim_LpSolver, regression_multiple_solutions){
cv::Mat A,B,z,etalon_z;
#if 1
//trivial example with multiple solutions
A=(cv::Mat_<double>(2,1)<<1,1);
B=(cv::Mat_<double>(1,3)<<1,1,1);
std::cout<<"here A goes\n"<<A<<"\n";
int res=cv::solveLP(A,B,z);
printf("res=%d\n",res);
printf("scalar %g\n",z.dot(A));
std::cout<<"here z goes\n"<<z<<"\n";
ASSERT_EQ(res,1);
ASSERT_EQ(z.dot(A),1);
#endif
}
TEST(Optim_LpSolver, regression_cycling){
cv::Mat A,B,z,etalon_z;
#if 1
//example with cycling from http://people.orie.cornell.edu/miketodd/or630/SimplexCyclingExample.pdf
A=(cv::Mat_<double>(4,1)<<10,-57,-9,-24);
B=(cv::Mat_<double>(3,5)<<0.5,-5.5,-2.5,9,0,0.5,-1.5,-0.5,1,0,1,0,0,0,1);
std::cout<<"here A goes\n"<<A<<"\n";
int res=cv::solveLP(A,B,z);
printf("res=%d\n",res);
printf("scalar %g\n",z.dot(A));
std::cout<<"here z goes\n"<<z<<"\n";
ASSERT_EQ(z.dot(A),1);
//ASSERT_EQ(res,1);
#endif
}

22
modules/core/test/test_misc.cpp
View File

@@ -3,28 +3,6 @@
using namespace cv;
using namespace std;
TEST(Core_Drawing, _914)
{
const int rows = 256;
const int cols = 256;
Mat img(rows, cols, CV_8UC1, Scalar(255));
line(img, Point(0, 10), Point(255, 10), Scalar(0), 2, 4);
line(img, Point(-5, 20), Point(260, 20), Scalar(0), 2, 4);
line(img, Point(10, 0), Point(10, 255), Scalar(0), 2, 4);
double x0 = 0.0/pow(2.0, -2.0);
double x1 = 255.0/pow(2.0, -2.0);
double y = 30.5/pow(2.0, -2.0);
line(img, Point(int(x0), int(y)), Point(int(x1), int(y)), Scalar(0), 2, 4, 2);
int pixelsDrawn = rows*cols - countNonZero(img);
ASSERT_EQ( (3*rows + cols)*3 - 3*9, pixelsDrawn);
}
TEST(Core_OutputArrayCreate, _1997)
{
struct local {