/* rimage_matlab.c
 * Created: Apr 3, 2006
 * Author: Andrew Salamon
 *
 * Convert from matlab arrays to RImage<float>'s and back again.
 * Write and read Matlab images to/from a mat file.
 */

#include "rimage_matlab.h"

/** Write one Matlab array to a file.
 * The array will be written out with the name "array", but will be the
 *  only variable in the file so the name should never be needed.
 * This will be used to store input and output data for automated testing.
 *
 * See also readImageMat() for getting the array back out of the file.
 */
void writeImageMat( const mxArray *img, const char *fileName )
{
	int status = 0;
	MATFile *fileID = matOpen( fileName, "w" );

  if ( fileID == NULL ) {
#ifdef DEBUG
    printf("Error creating file %s\n", fileName );
#endif
    return;
  }

  status = matPutVariable( fileID, "array", img );
  if (status != 0) {
#ifdef DEBUG
      printf("%s :  Error using matPutVariable on line %d\n", __FILE__, __LINE__);
#endif
      return;
  }

  if( matClose(fileID) != 0) {
#ifdef DEBUG
    printf("Error closing file %s\n",fileName);
#endif
    return;
  }
}

/** Read one Matlab array from a file.
 * The first variable in a Matlab .mat file will be read in and returned.
 * This is primarily used to create an RImage object from a saved .mat file and hasn't been tested for any other use.
 *
 * See also writeImageMat() for saving the array into the file.
 */
mxArray * readImageMat( const char *fileName )
{
	MATFile *fileID = matOpen( fileName, "r" );
	const char *name;
	mxArray *pa;

	if ( fileID == NULL ) {
#ifdef DEBUG
		printf("Error opening file %s for reading\n", fileName );
#endif
		return NULL;
	}

    pa = matGetNextVariable( fileID, &name);
    if( pa == NULL )
	{
		matClose( fileID );
		return NULL;
    }

	matClose(fileID);

	return pa;
}

/** Create an RImage from a Matlab array.
 * Takes an mxArray, possibly from readImageMat(), and returns an RImage.
 * The RImage is created with new and needs to be deleted when the caller is done with it.
 */
RImage<float> *rimageFromMatlabArray( const mxArray *array )
{
	double *img = mxGetPr(array);
	int rows = static_cast<int>(mxGetM(array));
	int cols = static_cast<int>(mxGetN(array));

	if( !rows || !cols )
		return NULL;

	// create image and make upright (i.e., transpose)
	RImage<float> *pixels = new RImage<float>(cols, rows);
	for (int c = 0; c < cols; c++)
	{
		for(int r = 0; r < rows; r++)
		{
			pixels->array[cols*r+c] = static_cast<float>(*(img++));
		}
	}

	return pixels;
}

/** Create a Matlab array from an RImage.
 * The mxArray is created with mxCreateDoubleMatrix and needs to be deleted when
 * the caller is done with it, unless it's being returned to Matlab from a mex.
 */
mxArray *matlabArrayFromRimage( const RImage<float> *image )
{
	mxArray *result;

    if( image->numpixels  == 0 )
	{
        result = mxCreateDoubleMatrix(0,0,mxREAL);
    } else {
		int cols = image->width;
		int rows = image->height;
        result = mxCreateDoubleMatrix( rows, cols, mxREAL );
        double *plx = mxGetPr( result );
        float *resPtr = image->array;

		for (int c = 0; c < cols; c++)
		{
			for(int r = 0; r < rows; r++)
			{
				plx[cols*r+c] = *resPtr++;
			}
		}
	}

	return result;
}

/** Create a Matlab array from an array of doubles.
 * The mxArray is created with mxCreateDoubleMatrix and needs to be deleted when
 * the caller is done with it, unless it's being returned to Matlab from a mex.
 */
mxArray *matlabArrayFromDoubles( const double *image, int width, int height )
{
	mxArray *result;
    int numpixels = width * height;

    if( numpixels  == 0 )
	{
        result = mxCreateDoubleMatrix(0,0,mxREAL);
    } else {
		int cols = width;
		int rows = height;
        result = mxCreateDoubleMatrix( rows, cols, mxREAL );
        double *plx = mxGetPr( result );
        const double *resPtr = image;

		for (int c = 0; c < cols; c++)
		{
			for(int r = 0; r < rows; r++)
			{
				plx[cols*r+c] = *resPtr++;
			}
		}
	}

	return result;
}

/** Create a Matlab array from an array of complex.
 * The mxArray is created with mxCreateDoubleMatrix and needs to be deleted when
 * the caller is done with it, unless it's being returned to Matlab from a mex.
 */
mxArray *matlabArrayFromComplex( const complex<double> *image, int width, int height )
{
	mxArray *result;
    int numpixels = width * height;

    if( numpixels  == 0 )
	{
        result = mxCreateDoubleMatrix(0,0,mxREAL);
    } else {
        result = mxCreateDoubleMatrix( height, width, mxCOMPLEX );
        double *plx = mxGetPr( result );
        double *pli = mxGetPi( result );
        const complex<double> *resPtr = image;

		for (int c = 0; c < width; c++)
		{
			for(int r = 0; r < height; r++)
			{
				plx[height*c+r] = resPtr->real(); // (*resPtr)[0];
				pli[height*c+r] = resPtr->imag(); // (*resPtr)[1];
                ++resPtr;
			}
		}
	}

	return result;
}

/** Create a Matlab array from an array of complex.
 * The mxArray is created with mxCreateDoubleMatrix and needs to be deleted when
 * the caller is done with it, unless it's being returned to Matlab from a mex.
 */
complex<double> *complexFromMatlabArray( const mxArray *array )
{
    double *real = mxGetPr(array);
    double *imag = mxGetPi(array);
    int rows = static_cast<int>(mxGetM(array));
    int cols = static_cast<int>(mxGetN(array));
    int numpixels = rows * cols;

    if( numpixels <= 0 )
	{
        return NULL;
    }
    else
    {
        complex<double> *res = (complex<double> *)fftw_malloc( sizeof(complex<double>) * numpixels );
        for(int c = 0; c < cols; c++)
        {
            for(int r = 0; r < rows; r++)
            {
                res[cols*r+c] = complex<double>( *(real++), *(imag++) );
            }
        }
        return res;
	}

	return NULL;
}

/** Create a Matlab array from an array of fftw complexes.
 * The mxArray is created with mxCreateDoubleMatrix and needs to be deleted when
 * the caller is done with it, unless it's being returned to Matlab from a mex.
 */
mxArray *matlabArrayFromComplexFFTW( const fftw_complex *image, int width, int height )
{
	mxArray *result;
    int numpixels = width * height;

    if( numpixels  == 0 )
	{
        result = mxCreateDoubleMatrix(0,0,mxREAL);
    } else {
		int num = width*height;
        result = mxCreateDoubleMatrix( height, width, mxCOMPLEX );
        double *plx = mxGetPr( result );
        double *pli = mxGetPi( result );
        const fftw_complex *resPtr = image;

		for (int i = 0; i < num; ++i)
        {
            plx[i] = (*resPtr)[0];
            pli[i] = (*resPtr)[1];
            ++resPtr;
        }
	}

	return result;
}