www.pudn.com > Image_segment.rar > Core.inl
//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
// * No commercial use is allowed.
// This software can only be used
// for non-commercial purposes. This
// distribution is mainly intended for
// academic research and teaching.
// * Redistributions of source code must
// retain the above copyright notice, this
// list of conditions and the following
// disclaimer.
// * Redistributions of binary form must
// mention the above copyright notice, this
// list of conditions and the following
// disclaimer in a clearly visible part
// in associated product manual,
// readme, and web site of the redistributed
// software.
// * Redistributions 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 Baris Sumengen 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
//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.
namespace MathCore
{
template< class T >
void Plot(Vector& v, int height)
{
int bufferWidth = 80;
Vector ind(bufferWidth-2);
for(int i=0; i v.Numel()-1)
{
newInd = v.Numel()-1;
}
ind(i) = v(newInd);
}
T max = Maximum(ind);
T min = Minimum(ind);
ind -= min;
ind *= (T)((height-1.1)/(max-min));
Matrix m(height, ind.Numel(), 0);
for(int i=0; i
Vector Find(Vector& m)
{
std::vector ind;
for(int i=0; i temp((int)ind.size());
memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
//std::vector::const_iterator constIterator;
//int j = 0;
//for(constIterator = ind.begin(); constIterator != ind.end(); constIterator++)
//{
// temp[j] = *constIterator;
// j++;
//}
return temp;
}
template
Vector Find(Matrix& m)
{
std::vector ind;
for(int i=0; i temp((int)ind.size());
memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
return temp;
}
template
void Find(Matrix& m, Vector& I, Vector& J)
{
std::vector indI;
std::vector indJ;
for(int j=0; j((int)indI.size());
memcpy(I.Data(), &indI[0], sizeof(int)*indI.size() );
J = Vector((int)indJ.size());
memcpy(J.Data(), &indJ[0], sizeof(int)*indJ.size() );
}
// Vectors need to be of length 3
template
Vector Cross(Vector& x, Vector& y)
{
Vector temp(3);
temp[0] = x[1]*y[2]-x[2]*y[1];
temp[1] = x[2]*y[0]-x[0]*y[2];
temp[2] = x[0]*y[1]-x[1]*y[0];
return temp;
}
template
Matrix Cross(Matrix& x, Matrix& y)
{
if(x.Rows() == 3)
{
return Cross(x,y,1);
}
else if(x.Columns() == 3)
{
return Cross(x,y,2);
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("One of the dimension of the matrices x and y should be of length 3!");
}
}
template
Matrix Cross(Matrix& x, Matrix& y, int dimension)
{
if(x.Rows() != y.Rows() || x.Columns() != y.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrices are not of the same size!");
}
Matrix temp(x.Rows(), x.Columns());
if(dimension == 1)
{
if(temp.Rows() != 3)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Both matrices x and y should be of length 3 along the dimension on which the cross product is taken!");
}
for(int z=0; z
Vector CumSum(Vector& m)
{
Vector temp(m.Length());
temp.ElemNC(0) = m.ElemNC(0);
for(int i=1; i
Matrix CumSum(Matrix& m)
{
return CumSum(m,1);
}
template
Matrix CumSum(Matrix& m, int dimension)
{
Matrix temp(m.Rows(), m.Columns());
if(dimension == 1)
{
for(int z=0; z
Vector& CumSumI(Vector& m)
{
for(int i=1; i
Matrix& CumSumI(Matrix& m)
{
return CumSum(m,1);
}
template
Matrix& CumSumI(Matrix& m, int dimension)
{
if(dimension == 1)
{
for(int z=0; z
Vector CumProd(Vector& m)
{
Vector temp(m.Length());
temp.ElemNC(0) = m.ElemNC(0);
for(int i=1; i
Matrix CumProd(Matrix& m)
{
return CumProd(m,1);
}
template
Matrix CumProd(Matrix& m, int dimension)
{
Matrix temp(m.Rows(), m.Columns());
if(dimension == 1)
{
for(int z=0; z
Vector& CumProdI(Vector& m)
{
for(int i=1; i
Matrix& CumProdI(Matrix& m)
{
return CumProd(m,1);
}
template
Matrix& CumProdI(Matrix& m, int dimension)
{
if(dimension == 1)
{
for(int z=0; z
Matrix Diag(Vector& m)
{
return Diag(m,0);
}
template
Matrix Diag(Vector& m, int k)
{
int absk = abs(k);
int dim = m.Length() + absk;
Matrix temp(dim,dim,0);
if(k>=0)
{
for(int i=0; i
Vector Diag(Matrix& m)
{
return Diag(m,0);
}
template
Vector Diag(Matrix& m, int k)
{
Vector temp;
int absk = abs(k);
int dim;
if(k>=0)
{
if(absk>=m.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
}
dim = m.Columns()-absk > m.Rows() ? m.Rows() : m.Columns()-absk;
temp = Vector(dim);
for(int i=0; i=m.Rows())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
}
dim = m.Rows()-absk > m.Columns() ? m.Columns() : m.Rows()-absk;
temp = Vector(dim);
for(int i=0; i
Matrix FlipLR(Matrix& m)
{
Matrix temp(m.Rows(), m.Columns());
for(int j=0; j
Matrix& FlipLRI(Matrix& m)
{
for(int j=0; j
Matrix FlipUD(Matrix& m)
{
Matrix temp(m.Rows(), m.Columns());
for(int j=0; j
Matrix& FlipUDI(Matrix& m)
{
for(int j=0; j
Vector Reverse(Vector& m)
{
Vector temp(m.Length());
for(int i=0; i
Vector& ReverseI(Vector& m)
{
for(int i=0; i<(m.Length()+1)/2; i++)
{
T dummy = m.ElemNC(i);
m.ElemNC(i) = m.ElemNC(m.Length()-i-1);
m.ElemNC(m.Length()-i-1) = dummy;
}
return m;
}
template
Matrix Rot90(Matrix& m)
{
return Rot90(m,1);
}
template
Matrix Rot90(Matrix& m, int k)
{
int rot = k % 4;
if(rot<0)
{
rot += 4;
}
if(rot==0)
{
return m.Clone();
}
else if(rot==1)
{
return FlipUD(~m);
}
else if(rot==2)
{
return FlipLR(FlipUD(m));
}
else if(rot==3)
{
return FlipLR(~m);
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Something went wrong when trying to rotate the matrix!");
}
}
template
Vector Maximum(Matrix& m)
{
Vector maxes(m.Columns());
for(int j=0; j
T Maximum(Vector& m)
{
T max = m.ElemNC(0);
for(int i=1; i= m.ElemNC(i) ? max : m.ElemNC(i);
}
return max;
}
template
Vector Minimum(Matrix& m)
{
Vector mins(m.Columns());
for(int j=0; j
T Minimum(Vector& m)
{
T min = m.ElemNC(0);
for(int i=1; i
Vector Mean(Matrix& m)
{
Vector means(m.Columns());
for(int j=0; j
double Mean(Vector& m)
{
double mean = 0;
for(int i=0; i
Vector Median(Matrix& m)
{
Vector medians(m.Columns());
for(int j=0; j
T Median(Vector& m)
{
Vector temp = m.Clone();
T *arr = temp.Data();
int n = temp.Length();
int low, high;
int median;
int middle, ll, hh;
low = 0;
high = m.Length()-1;
median = (low + high)/2;
for(;;)
{
if(high <= low)
{
return arr[median];
}
if(high == low + 1)
{
if(arr[low] > arr[high])
{
Swap(arr[low], arr[high]);
}
return arr[median];
}
middle = (low + high)/2;
if(arr[middle] > arr[high])
{
Swap(arr[middle], arr[high]);
}
if(arr[low] > arr[high])
{
Swap(arr[low], arr[high]);
}
if(arr[middle] > arr[low])
{
Swap(arr[middle], arr[low]);
}
Swap(arr[middle], arr[low+1]) ;
ll = low + 1;
hh = high;
for(;;)
{
do ll++; while(arr[low] > arr[ll]);
do hh--; while(arr[hh] > arr[low]);
if (hh < ll)
{
break;
}
Swap(arr[ll], arr[hh]);
}
Swap(arr[low], arr[hh]);
if (hh <= median)
{
low = ll;
}
if (hh >= median)
{
high = hh - 1;
}
}
}
template
Vector Var(Matrix& m)
{
return Var(m,true);
}
template
double Var(Vector& m)
{
return Var(m,true);
}
template
Vector Var(Matrix& m, bool unbiased)
{
Vector vars(m.Columns());
for(int j=0; j
double Var(Vector& m, bool unbiased)
{
double mean = Mean(m);
double var = 0;
for(int i=0; i
Vector Std(Matrix& m)
{
return Std(m,true);
}
template
double Std(Vector& m)
{
return Std(m,true);
}
template
Vector Std(Matrix& m, bool unbiased)
{
Vector stds(m.Columns());
for(int j=0; j
double Std(Vector& m, bool unbiased)
{
return sqrt(Var(m,unbiased));
}
template
Vector Sum(Matrix& m)
{
return Sum(m,1);
}
template
Vector Sum(Matrix& m, int dimension)
{
Vector sums;
if(dimension == 1)
{
sums = Vector(m.Columns());
for(int j=0; j(m.Rows());
for(int j=0; j
T Sum(Vector& m)
{
T sum = 0;
for(int i=0; i
Vector Prod(Matrix& m)
{
return Prod(m,1);
}
template
Vector Prod(Matrix& m, int dimension)
{
Vector prods;
if(dimension == 1)
{
prods = Vector(m.Columns());
for(int j=0; j(m.Rows());
for(int j=0; j
T Prod(Vector& m)
{
T prod = 1;
for(int i=0; i
Vector Diff(Vector& m)
{
if(m.Length() == 1)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Diff can only be called on vectors of length 2 or more!");
}
Vector temp(m.Length()-1);
for(int i=0;i
Matrix Diff(Matrix& m)
{
if(m.Rows() == 1)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Diff can only be called on matrices with 2 or more rows!");
}
Matrix temp(m.Rows()-1, m.Columns());
for(int j=0; j
Matrix RepMat(Matrix& m, int M, int N)
{
Matrix temp(M*m.Rows(),N*m.Columns());
for(int k=0; k
Matrix Reshape(Matrix& m, int M, int N)
{
if(m.Length() != M*N)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
}
Matrix temp(M,N);
memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
return temp;
}
template
Matrix Reshape(Vector& m, int M, int N)
{
if(m.Length() != M*N)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
}
Matrix temp(M,N);
memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
return temp;
}
template
Vector& QuickSort(Vector& m)
{
int M = 7;
int NSTACK = 50;
int i, j, k;
int l=1;
int ir = m.Length();
Vector stack(m.Length());
int *istack = stack.Data()-1;
int jstack = 0;
T a;
T* arr = m.Data()-1;
for (;;)
{
if (ir-l < M)
{
for (j=l+1;j<=ir;j++)
{
a=arr[j];
for (i=j-1;i>=l;i--)
{
if (arr[i] <= a)
{
break;
}
arr[i+1]=arr[i];
}
arr[i+1]=a;
}
if (jstack == 0)
{
break;
}
ir=istack[jstack--];
l=istack[jstack--];
}
else
{
k=(l+ir) >> 1;
Swap(arr[k],arr[l+1]);
if (arr[l] > arr[ir])
{
Swap(arr[l],arr[ir]);
}
if (arr[l+1] > arr[ir])
{
Swap(arr[l+1],arr[ir]);
}
if (arr[l] > arr[l+1])
{
Swap(arr[l],arr[l+1]);
}
i=l+1;
j=ir;
a=arr[l+1];
for (;;)
{
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i)
{
break;
}
Swap(arr[i],arr[j]);
}
arr[l+1]=arr[j];
arr[j]=a;
jstack += 2;
if (jstack > NSTACK)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Stack size is too small in quick sort!");
}
if (ir-i+1 >= j-l)
{
istack[jstack]=ir;
istack[jstack-1]=i;
ir=j-1;
}
else
{
istack[jstack]=j-1;
istack[jstack-1]=l;
l=i;
}
}
}
return m;
}
template
Vector Sort(Vector& m)
{
Vector sorted = m.Clone();
QuickSort(sorted);
return sorted;
}
template
Matrix Sort(Matrix& m, int dimension)
{
Matrix temp;
if(dimension == 1)
{
temp = m.Clone();
for(int j=0; j(m.Rows(), m.Columns());
for(int i=0; i dummy(temp.Columns());
for(int j=0; j
Matrix Sort(Matrix& m)
{
return Sort(m,1);
}
template
Vector& SortI(Vector& m)
{
return QuickSort(m);
}
template
Matrix& SortI(Matrix& m, int dimension)
{
if(dimension == 1)
{
for(int j=0; j dummy(m.Columns());
for(int i=0; i
Matrix& SortI(Matrix& m)
{
return SortI(m,1);
}
template
Vector Abs(Vector& m)
{
Vector dummy(m.Numel());
for(int i=0; i
Vector& AbsI(Vector& m)
{
for(int i=0; i
Matrix Abs(Matrix& m)
{
Matrix dummy(m.Rows(), m.Columns());
for(int i=0; i
Matrix& AbsI(Matrix& m)
{
for(int i=0; i
void MeshGrid(Vector