diff --git a/modules/core/src/matrix.cpp b/modules/core/src/matrix.cpp
index aaaac10..fcd7ba3 100644
--- a/modules/core/src/matrix.cpp
+++ b/modules/core/src/matrix.cpp
@@ -2399,6 +2399,49 @@ static void generateRandomCenter(const vector<Vec2f>& box, float* center, RNG& r
         center[j] = ((float)rng*(1.f+margin*2.f)-margin)*(box[j][1] - box[j][0]) + box[j][0];
 }
 
+#ifdef HAVE_TBB 
+
+class KMeansPPDistanceComputer
+{
+  float *tdist2;
+  const float *data;
+  const float *dist;
+  const int dims;
+  const int step;
+  const int stepci;
+
+public:
+
+  KMeansPPDistanceComputer(float *_tdist2,
+			   const float *_data,
+			   const float *_dist,
+			   int _dims,
+			   int _step,
+			   int _stepci)
+    : tdist2(_tdist2),
+      data(_data),
+      dist(_dist),
+      dims(_dims),
+      step(_step),
+      stepci(_stepci)
+  {
+
+  }
+
+  void operator()(const cv::BlockedRange& range) const
+  {
+
+    const int begin = range.begin();
+    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]);
+    }
+  }
+  
+};
+
+#endif //HAVE_TBB
 
 /*
 k-means center initialization using the following algorithm:
@@ -2436,11 +2479,21 @@ static void generateCentersPP(const Mat& _data, Mat& _out_centers,
                 if( (p -= dist[i]) <= 0 )
                     break;
             int ci = i;
+#ifndef HAVE_TBB
             for( i = 0; i < N; i++ )
             {
                 tdist2[i] = std::min(normL2Sqr_(data + step*i, data + step*ci, dims), dist[i]);
                 s += tdist2[i];
             }
+#else
+	    cv::parallel_for(cv::BlockedRange(0, N),
+			     KMeansPPDistanceComputer(tdist2, data, dist, dims, step, step*ci));
+	    
+            for( i = 0; i < N; i++ )
+            {
+	      s += tdist2[i];
+            }
+#endif //HAVE_TBB
 
             if( s < bestSum )
             {
@@ -2463,6 +2516,67 @@ static void generateCentersPP(const Mat& _data, Mat& _out_centers,
     }
 }
 
+    
+#ifdef HAVE_TBB 
+
+class KMeansDistanceComputer
+{
+  float *distances;
+  int *labels;
+  const cv::Mat& data;
+  const cv::Mat& centers;
+  
+public:
+
+  KMeansDistanceComputer(float *_distances,
+			 int *_labels,
+			 const cv::Mat& _data,
+			 const cv::Mat& _centers)
+    : distances(_distances),
+      labels(_labels),
+      data(_data),
+      centers(_centers)
+  {
+    CV_DbgAssert(centers.cols == data.cols);
+  }
+
+  void operator()(const cv::BlockedRange& range) const
+  {
+    const int begin = range.begin();
+    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;
+      }
+    
+  }
+  
+};
+
+#endif //HAVE_TBB
+
 }
 
 double cv::kmeans( InputArray _data, int K,
@@ -2508,6 +2622,10 @@ double cv::kmeans( InputArray _data, int K,
     vector<Vec2f> _box(dims);
     Vec2f* box = &_box[0];
 
+#ifdef HAVE_TBB 
+    cv::Mat dists(N, 1, CV_32F);
+#endif //HAVE_TBB
+
     double best_compactness = DBL_MAX, compactness = 0;
     RNG& rng = theRNG();
     int a, iter, i, j, k;
@@ -2682,6 +2800,7 @@ double cv::kmeans( InputArray _data, int K,
                 break;
 
             // assign labels
+#ifndef HAVE_TBB 
             compactness = 0;
             for( i = 0; i < N; i++ )
             {
@@ -2704,6 +2823,17 @@ double cv::kmeans( InputArray _data, int K,
                 compactness += min_dist;
                 labels[i] = k_best;
             }
+#else
+	    float* dist = dists.ptr<float>(0);
+	    cv::parallel_for(cv::BlockedRange(0, N),
+			     KMeansDistanceComputer(dist, labels, data, centers));
+	    compactness = 0;
+	    for( i = 0; i < N; i++ )
+	    {
+	      compactness += dist[i];
+	    }
+
+#endif //HAVE_TBB
         }
 
         if( compactness < best_compactness )
