surfarray.cpp

A set of C++ and Matlab routines implementing the surfacelet transform surfacelet的一个非常好用的工具箱

		nRows *= pDims_logical[i];
	}

	double *pMagnitude = pData, *pMagMirror;
	double *pReal = pData, *pImag = pData + 1;
	lengthRow = pDims_real[nDims - 1] / 2;

	for (i = 0; i < nRows; i++)
	{
		for (j = 0; j < lengthRow; j++)
		{
			*(pMagnitude++) = sqrt((*pReal) * (*pReal) + (*pImag) * (*pImag));
			pReal += 2;
			pImag += 2;
		}
	
		// tricky!
		// verify the following by trying two cases
		// (1) pDims_logicalp[nDims - 1] is even
		// (2) pDims_logicalp[nDims - 1] is odd
		pMagMirror = pMagnitude - padding;
		for (j = 0; j < lengthRow - padding; j++)
			*(pMagnitude++) = *(pMagMirror--);

		pMagnitude += padding;
		
	}

	IsReal = true;
}


//////////////////////////////////////////////////////////////////////////
//  Fill the current array with zeros
//////////////////////////////////////////////////////////////////////////

void SurfArray::ZeroFill()
{
    if (pData)
    {
        int i, n_total_numbers = 1;
        for (i = 0; i < nDims; i++)
            n_total_numbers *= pDims_real[i];

        memset(pData, 0, n_total_numbers * sizeof(double));
    }
}

//////////////////////////////////////////////////////////////////////////
//	Repeat the current array along each of its dimensions
//////////////////////////////////////////////////////////////////////////
//
//	WARNING:
//		
//	This routine has not be tested.

void RepeatArray(int N, int pDims[], int repRatio[], double *src, double *dst, int padding, double *pBuffer)
{
    double *pVal, *pVal_mirror;
    register int i;
    int nMirror;
    if (N == 1)
    {
        if (repRatio[0] == 1)
            memcpy(dst, src, (pDims[0] + padding) * sizeof(double));
        else
        {
            memcpy(pBuffer, src, (pDims[0] + padding) * sizeof(double));
            pVal = pBuffer + pDims[0] + padding;
            pVal_mirror = pVal - padding * 2;
            nMirror = (pDims[0] - padding) / 2;
            for (i = 0; i < nMirror; i++)
            {
                *(pVal++) = *(pVal_mirror++);
                *(pVal++) = -*pVal_mirror;
                pVal_mirror -= 2;
            }

            int targetLen = pDims[0] * repRatio[0];
            if (targetLen % 2)
                targetLen += 1;
            else
                targetLen += 2;

            int bufferLen = pDims[0] * 2;
            while (targetLen > 0)
            {
                memcpy(dst, pBuffer, bufferLen * sizeof(double));
                dst += bufferLen;
                targetLen -= bufferLen;
            }
        }
    }
    else
    {
        int shift_src = 1, shift_dst = 1;
        for (i = 1; i < N; i++)
        {
            shift_src *= pDims[i];
            shift_dst *= (pDims[i] * repRatio[i]);
        }

        for (i = 0; i < pDims[0]; i++)
        {
            RepeatArray(N - 1, pDims + 1, repRatio + 1, src + shift_src * i, dst + shift_dst * i, padding, pBuffer);
        }

        shift_dst *= pDims[0];

        for (i = 0; i < repRatio[0]; i++)
            memcpy(dst + shift_dst * (i + 1), dst, shift_dst * sizeof(double));
    }
}


//////////////////////////////////////////////////////////////////////////
//	Upsample the current array (in the frequency domain)
//////////////////////////////////////////////////////////////////////////
// 
//	PARAMETERS:
//
//	newArray
//		pointer to the upsampled array
// 
//	upRatio
//		upsampling ratio along each dimension
//
//	WARNING:
//		This routine has not been tested.

void SurfArray::UpsampleF(SurfArray* newArray, int upRatio[])
{
    // The original array must be in the frequency domain
    assert(!IsReal);

    int *pDimsNew;
    double *pBuffer;
    try
    {
        pDimsNew = new int [nDims];
        pBuffer = new double [pDims_logical[nDims - 1] * 2];
    }
    catch (std::bad_alloc)
    {
    	MessageLog("SurfArray", "UpsampleF", "Out of memory!");
        if (pDimsNew) delete [] pDimsNew;
        throw;
    }

    newArray->GetDims(pDimsNew);
    int i, checkDims = 0, padding;
    for (i = 0; i < nDims; i++)
        checkDims += abs(pDims_logical[i] * upRatio[i] - pDimsNew[i]);

    assert(checkDims == 0);

    // newArray must be in the Fourier domain
    newArray->SetIsRealValued(false);

    double *dst;
    dst = newArray->GetPointer(padding);
    // Similar to "repmat" in MATLAB.
    RepeatArray(nDims, pDims_logical, upRatio, pData, dst, padding, pBuffer);
    
    delete [] pBuffer;
    delete [] pDimsNew;
   
}


//////////////////////////////////////////////////////////////////////////
//	Multiply the current array with a 2-D slice along a certian signal plane
//////////////////////////////////////////////////////////////////////////

void SurfArray::Multiply2dComplexSlice(int dim1, int dim2, SurfMatrix* pSlice)
{
	// Not implemented yet.
	// This function is not used in the current Surfacelet implementation
	assert(false);
}


//////////////////////////////////////////////////////////////////////////
//	Export the current array (in the frequency domain) with Hermitian
//	symmetry along a specified dimension
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//
//	NewAxis
//		The new axis along which the Hermitian symmetry is utilized

double* SurfArray::NewHermitianSymmetryAxis(int NewAxis)
{
    // we are currently utilizing the symmetry along the last dimension
	// so the NewAxis must be different from that.
    assert((NewAxis <= nDims - 1) && (NewAxis >= 0));
    assert(!IsReal);
	assert(pDims_logical[NewAxis] > 2);
    
    double *pOut;
    int *pNewDims, *pIndices, *pMemoryDistance;
    int n_total_points = 1;
    register int n, m, j, i;

    try
    {
        pNewDims = new int [nDims];
        pIndices = new int [nDims];
        pMemoryDistance = new int [nDims-1];
        for (i = 0; i < nDims; i++)
            pNewDims[i] = pDims_logical[i];
        pNewDims[NewAxis] = pNewDims[NewAxis] / 2 + 1;

        for (i = 0; i < nDims; i++)
            n_total_points *= pNewDims[i];
        // 1 complex = 2 double
		n_total_points *= 2;        

        pOut = new double [n_total_points];
        
    }
    catch (std::bad_alloc)
    {
    	MessageLog("SurfArray", "NewHermitianSymmetryAxis", "Out of memory!");
        throw;
    }

    if (NewAxis == nDims - 1)
    {
        memcpy(pOut, pData, n_total_points * sizeof(double));
        delete [] pMemoryDistance;    
        delete [] pIndices;
        delete [] pNewDims;
        // IMPORTANT
        return pOut;
    }

    int dimA = 1;
    for (n = 0; n < NewAxis; n++)
        dimA *= pDims_logical[n];

    int dimB = pNewDims[NewAxis];

    int dimC = 1;
    for (n = NewAxis + 1; n < nDims - 1; n++)
        dimC *= pDims_logical[n];

    int dimD = pDims_logical[nDims - 1] - (pDims_logical[nDims - 1] / 2 + 1);

	// Initialize the indices
    for (n = 0; n < nDims; n++)
        pIndices[n] = 0;
    
    pMemoryDistance[nDims - 2] = pDims_real[nDims-1];
    for (n = nDims - 3; n >= 0; n--)
        pMemoryDistance[n] = pMemoryDistance[n+1] * pDims_real[n+1];


    double *dst, *src, *src_symm;
    dst = pOut;
    src = pData;

    int skip = pDims_real[nDims - 1];
    int skip_symm, offset;
	int src_skip = pMemoryDistance[NewAxis] * (pDims_logical[NewAxis]-dimB);
	offset = 2 * pDims_logical[nDims - 1] - pDims_real[nDims - 1];

    for (i = 0; i < dimA; i++)
	{
        for (j = 0; j < dimB; j++)
            for (m = 0; m < dimC; m++)
            {
                memcpy(dst, src, skip * sizeof(double));
                dst += skip;
                src += skip;
                
                // update src_symm
                skip_symm = 0;
                for (n = 0; n < nDims - 1; n++)
					if (pIndices[n])
						skip_symm += (pDims_logical[n] - pIndices[n]) * pMemoryDistance[n];
                
				skip_symm += offset;

                src_symm = pData + skip_symm;

                for (n = 0; n < dimD; n++)
                {
                    *(dst++) = *(src_symm++);
                    *(dst++) = -*(src_symm++);
                    src_symm -= 4;
                }

                // update indices
                // Update pIndices
                for (n = nDims - 2; n >= 0; n--)
                {
                    if ( (++pIndices[n]) < pNewDims[n])
                        break;
                    else
                        pIndices[n] = 0;

                }
            }
		src += src_skip;
	}
        
    delete [] pMemoryDistance;    
    delete [] pIndices;
    delete [] pNewDims;

    return pOut;
}


//////////////////////////////////////////////////////////////////////////
//	Import frequency data into the current array
//////////////////////////////////////////////////////////////////////////
//
//	PARAMETERS:
//	
//	pIn
//		The pointer to the outsize memory space
//
//	Previous
//		The axis along which the original data utilize the Hermitian symmetry

void SurfArray::RestoreHermianSymmetryAxis(double *pIn, int PreviousAxis)
{
	// we are currently utilizing the symmetry along the last dimension
	assert((PreviousAxis <= nDims - 1) && (PreviousAxis >= 0));
	assert(!IsReal);

    int *pPrevioisDims, *pIndices, *pMemoryDistance;
	int n_total_points = 1;
	register int n, m, j, i;

    int nTotalPoints = 1;
    if (PreviousAxis == nDims - 1)
    {
        for (n = 0; n < nDims; n++)
            nTotalPoints *= pDims_real[n];

        memcpy(pData, pIn, nTotalPoints * sizeof(double));

        // IMPORTANT
        return;
    }
   
	try
	{
		pPrevioisDims = new int [nDims];
		pIndices = new int [nDims];
		pMemoryDistance = new int [nDims-1];
		for (i = 0; i < nDims; i++)
			pPrevioisDims[i] = pDims_logical[i];
		pPrevioisDims[PreviousAxis] = pPrevioisDims[PreviousAxis] / 2 + 1;
	}
	catch (std::bad_alloc)
	{
		MessageLog("SurfArray", "NewHermitianSymmetryAxis", "Out of memory!");
		throw;
	}

	int dimA = 1;
	for (n = 0; n < PreviousAxis; n++)
		dimA *= pDims_logical[n];

	int dimB = pPrevioisDims[PreviousAxis];
	int dimB2 = pDims_logical[PreviousAxis] - dimB;

	int dimC = 1;
	for (n = PreviousAxis + 1; n < nDims - 1; n++)
		dimC *= pDims_logical[n];

	int dimD = pDims_real[nDims - 1] / 2;

	// Initialize the indices
	for (n = 0; n < nDims; n++)
		pIndices[n] = 0;

	pMemoryDistance[nDims - 2] = pPrevioisDims[nDims-1] * 2;
	for (n = nDims - 3; n >= 0; n--)
		pMemoryDistance[n] = pMemoryDistance[n+1] * pPrevioisDims[n+1];


	double *dst, *src, *src_symm;
	dst = pData;
	src = pIn;

	int skip = pDims_real[nDims - 1];
	int skip_symm, offset;
	offset = 2 * pDims_logical[nDims - 1] - 2;

	for (i = 0; i < dimA; i++)
	{
		for (j = 0; j < dimB; j++)
		{
			for (m = 0; m < dimC; m++)
			{
				memcpy(dst, src, skip*sizeof(double));
				dst += skip;
				src += pMemoryDistance[nDims - 2];
			}
		}

		pIndices[PreviousAxis] = dimB;

		for (j = 0; j < dimB2; j++)
		{
			for (m = 0; m < dimC; m++)
			{
				
				// update src_symm
				skip_symm = 0;
				for (n = 0; n < nDims - 1; n++)
					if (pIndices[n])
						skip_symm += (pDims_logical[n] - pIndices[n]) * pMemoryDistance[n];

				src_symm = pIn + skip_symm;
				*(dst++) = *src_symm;
				*(dst++) = -*(src_symm+ 1);

				src_symm += offset;

				for (n = 1; n < dimD; n++)
				{
					*(dst++) = *(src_symm++);
					*(dst++) = -*(src_symm++);
					src_symm -= 4;
				}

				// Update pIndices
				for (n = nDims - 2; n >= 0; n--)
				{
					if ( (++pIndices[n]) < pDims_logical[n])
						break;
					else
						pIndices[n] = 0;

				}
			}
		}
	}

	delete [] pMemoryDistance;    
	delete [] pIndices;
	delete [] pPrevioisDims;
}


//	This software is provided "as-is", without any express or implied
//	warranty. In no event will the authors be held liable for any 
//	damages arising from the use of this software.