www.pudn.com > Image_segment.rar > Analysis.cpp
//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.
#include "./Analysis.h"
#include "mkl.h"
namespace Analysis
{
void StatusDisplay(long status)
{
long classError;
classError = DftiErrorClass(status, DFTI_ERROR_CLASS);
if(! classError)
{
cerr << "Error Status is not a member of Predefined Error Class\n";
}
else
{
char* errorMessage = DftiErrorMessage(status);
cerr << "Error_message = " << errorMessage << endl;
}
}
Vector FFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Vector& FFTI(Vector& x)
{
int n = x.Length();
//Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Vector FFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Vector& FFTI(Vector& x)
{
int n = x.Length();
//Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Vector FFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_REAL, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
for(int i=1;i<(n+1)/2;i++)
{
y[n-i] = conj(y[i]);
}
return y;
}
Vector FFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_REAL, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
for(int i=1;i<(n+1)/2;i++)
{
y[n-i] = conj(y[i]);
}
return y;
}
// ////////////////////
// Backward transform
// ////////////////////
Vector IFFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Backward transform
status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Vector& IFFTI(Vector& x)
{
int n = x.Length();
//Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Backward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Vector IFFT(Vector& x)
{
int n = x.Length();
Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Backward transform
status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Vector& IFFTI(Vector& x)
{
int n = x.Length();
//Vector y(n);
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Backward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
// /////////////////
// 2D Transform
// /////////////////
Matrix FFT2(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Matrix& FFT2I(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
//Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Matrix FFT2(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Matrix& FFT2I(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
//Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeForward( DescHandle, x.Data() );
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Matrix IFFT2(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Matrix& IFFT2I(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
//Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Matrix IFFT2(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Matrix& IFFT2I(Matrix& x)
{
int rows = x.Rows();
int cols = x.Columns();
//Matrix y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
//=========================================================================
// Convolution
//=========================================================================
Vector Conv(Vector& input, Vector& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
}
Vector Conv(Vector& input, Vector& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
}
Vector Conv(Vector& input, Vector& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Length() > input.Length())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter length should not be larger than the input signal!");
}
int fftlength = input.Length()+filter.Length()-1;
if(PowerOfTwo)
{
fftlength = NextPow2(fftlength);
}
Vector temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Vector(fftlength);
for(int i=0; i c1 = Reverse( input.Slice(0, wid1-1) );
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
temp1.ReadFromVector(c2, input.Length());
}
}
else if(b == Replicate)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector c1(wid1, input.ElemNC(0));
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector c2(wid2, input.ElemNC(input.Length()-1));
temp1.ReadFromVector(c2, input.Length());
}
}
Vector temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Vector(fftlength);
for(int i=0; i(input.Length());
for(int i=0; i o = IFFT( FFT(temp1) * FFT(temp2) );
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Length()+filter.Length()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
Vector Conv(Vector& input, Vector& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Length() > input.Length())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter length should not be larger than the input signal!");
}
int fftlength = input.Length()+filter.Length()-1;
if(PowerOfTwo)
{
fftlength = NextPow2(fftlength);
}
Vector temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Vector(fftlength);
for(int i=0; i c1 = Reverse( input.Slice(0, wid1-1) );
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
temp1.ReadFromVector(c2, input.Length());
}
}
else if(b == Replicate)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector c1(wid1, input.ElemNC(0));
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector c2(wid2, input.ElemNC(input.Length()-1));
temp1.ReadFromVector(c2, input.Length());
}
}
Vector temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Vector(fftlength);
for(int i=0; i(input.Length());
for(int i=0; i o = IFFT( FFT(temp1) * FFT(temp2) );
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Length()+filter.Length()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
// =============
// Conv2
// =============
Matrix Conv2(Matrix& input, Matrix& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
}
Matrix Conv2(Vector& filter1, Vector& filter2, Matrix& input, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexFloat(filter1), ToComplexFloat(filter2), ToComplexFloat(input), fa, b, PowerOfTwo) );
}
Matrix Conv2(Matrix& input, Matrix& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
}
Matrix Conv2(Vector& filter1, Vector& filter2, Matrix& input, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexDouble(filter1), ToComplexDouble(filter2), ToComplexDouble(input), fa, b, PowerOfTwo) );
}
Matrix Conv2(Matrix& input, Matrix& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
}
// Allocate area for input using the border
int fftrows = input.Rows()+filter.Rows()-1;
int fftcols = input.Columns()+filter.Columns()-1;
if(PowerOfTwo)
{
fftrows = NextPow2(fftrows);
fftcols = NextPow2(fftcols);
}
Matrix temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Matrix(fftrows, fftcols,0);
for(int j=0; j c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
}
if(wid2 != 0)
{
Matrix c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
temp1.ReadFromMatrix(c2, 0, input.Columns());
}
int hei1 = filter.Rows()/2;
int hei2 = filter.Rows() - hei1 - 1;
if(hei1 != 0)
{
Matrix r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
}
if(hei2 != 0)
{
Matrix r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
temp1.ReadFromMatrix(r2, input.Rows(), 0);
}
if(hei1 != 0 && wid1 != 0)
{
Matrix q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
}
if(hei1 != 0 && wid2 != 0)
{
Matrix q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
}
if(hei2 != 0 && wid1 != 0)
{
Matrix q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
}
if(hei2 != 0 && wid2 != 0)
{
Matrix q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
}
}
else if(b == Replicate)
{
int wid1 = filter.Columns()/2;
int wid2 = filter.Columns() - wid1 - 1;
if(wid1 != 0)
{
Matrix c1 = input.Slice(0, input.Rows()-1, 0, 0);
temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
}
if(wid2 != 0)
{
Matrix c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
}
int hei1 = filter.Rows()/2;
int hei2 = filter.Rows() - hei1 - 1;
if(hei1 != 0)
{
Matrix r1 = input.Slice(0, 0, 0, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
}
if(hei2 != 0)
{
Matrix r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
}
if(hei1 != 0 && wid1 != 0)
{
Matrix q11( hei1, wid1, input.ElemNC(0,0) );
temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
}
if(hei1 != 0 && wid2 != 0)
{
Matrix q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
}
if(hei2 != 0 && wid1 != 0)
{
Matrix q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
}
if(hei2 != 0 && wid2 != 0)
{
Matrix q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
}
}
// Allocate area for the filter
Matrix temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Matrix(fftrows, fftcols,0);
for(int j=0; j(input.Rows(), input.Columns(),0);
for(int j=0; j o = IFFT2( FFT2(temp1)*FFT2(temp2) );
// Based on the filter area, crop the border.
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
Matrix Conv2(Vector& filter1, Vector& filter2, Matrix& input, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter sizes is smaller than the input size
if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
}
Vector *cols = new (std::nothrow) Vector[input.Columns()];
Utility::CheckPointer(cols);
for(int i=0; i *rows = new (std::nothrow) Vector[rlength];
Utility::CheckPointer(rows);
for(int i=0; i(input.Columns());
for(int j=0; j temp(rlength, clength);
for(int i=0; i Conv2(Matrix& input, Matrix& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
}
// Allocate area for input using the border
int fftrows = input.Rows()+filter.Rows()-1;
int fftcols = input.Columns()+filter.Columns()-1;
if(PowerOfTwo)
{
fftrows = NextPow2(fftrows);
fftcols = NextPow2(fftcols);
}
Matrix temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Matrix(fftrows, fftcols,0);
for(int j=0; j c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
}
if(wid2 != 0)
{
Matrix c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
temp1.ReadFromMatrix(c2, 0, input.Columns());
}
int hei1 = filter.Rows()/2;
int hei2 = filter.Rows() - hei1 - 1;
if(hei1 != 0)
{
Matrix r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
}
if(hei2 != 0)
{
Matrix r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
temp1.ReadFromMatrix(r2, input.Rows(), 0);
}
if(hei1 != 0 && wid1 != 0)
{
Matrix q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
}
if(hei1 != 0 && wid2 != 0)
{
Matrix q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
}
if(hei2 != 0 && wid1 != 0)
{
Matrix q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
}
if(hei2 != 0 && wid2 != 0)
{
Matrix q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
}
}
else if(b == Replicate)
{
int wid1 = filter.Columns()/2;
int wid2 = filter.Columns() - wid1 - 1;
if(wid1 != 0)
{
Matrix c1 = input.Slice(0, input.Rows()-1, 0, 0);
temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
}
if(wid2 != 0)
{
Matrix c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
}
int hei1 = filter.Rows()/2;
int hei2 = filter.Rows() - hei1 - 1;
if(hei1 != 0)
{
Matrix r1 = input.Slice(0, 0, 0, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
}
if(hei2 != 0)
{
Matrix r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
}
if(hei1 != 0 && wid1 != 0)
{
Matrix q11( hei1, wid1, input.ElemNC(0,0) );
temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
}
if(hei1 != 0 && wid2 != 0)
{
Matrix q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
}
if(hei2 != 0 && wid1 != 0)
{
Matrix q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
}
if(hei2 != 0 && wid2 != 0)
{
Matrix q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
}
}
// Allocate area for the filter
Matrix temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Matrix(fftrows, fftcols,0);
for(int j=0; j(input.Rows(), input.Columns(),0);
for(int j=0; j o = IFFT2( FFT2(temp1) * FFT2(temp2) );
// Based on the filter area, crop the border.
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
Matrix Conv2(Vector& filter1, Vector& filter2, Matrix& input, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter sizes is smaller than the input size
if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
}
Vector *cols = new (std::nothrow) Vector[input.Columns()];
Utility::CheckPointer(cols);
for(int i=0; i *rows = new (std::nothrow) Vector[rlength];
Utility::CheckPointer(rows);
for(int i=0; i(input.Columns());
for(int j=0; j temp(rlength, clength);
for(int i=0; i