matwrite.c

/*
  matread.c: Provides routines for outputting data to Level 5 .mat files.
  Copyright (C) 2005  Kevin McHale
  
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
*/
#include "matio_internal_types.h"
#include "matio.h"

int matio_write(const char * filename, MATdata * data)
{
	FILE *outfile;
	char buffer[116];
	unsigned char *compress_buffer;
	time_t t;
	MATdata * temp;
	int i, offset, status, *bigbuffer; 
	unsigned long compress_size, datasize;
	outfile = (FILE *) fopen(filename, "w");
	if(!outfile)
		perror(filename);

	strcpy(buffer, "MATLAB 5.0 MAT-file, Generated by MATIO (C) Kevin McHale 2005, Created on: ");
	time(&t);
	strcat(buffer, ctime(&t));

	memset(buffer + strlen(buffer), ' ', 116 - strlen(buffer));

	fwrite(buffer, 1, 116, outfile);

	memset(buffer, ' ', 8);

	fwrite(buffer, 1, 8, outfile);
	
	*((short int *) buffer) = 0x0100;
	fwrite(buffer, 2, 1, outfile);

	// Write the endian indicator:
	*((short int *) buffer) = 256 * (short int) 'M' + (short int) 'I';
	fwrite(buffer, 2, 1, outfile);

	temp = data;
	while (temp != NULL) {
		// Ignore any data types that we do not understand:	
		if(temp -> type != MATERROR) {
			// We first add up the entire size of the matrix, in bytes:		
			datasize = 0;
			// for the data subelement:
			datasize += (matio_array_size(temp) * 8 + 8) * ((temp -> type == MATCOMPLEX) ? 2 : 1); 
			// for the matrix header and array flags:
			datasize += 24; 
			// for the dimensions:
			datasize += 8 + temp -> num_dim * 4 + ((temp -> num_dim % 2) ? 4 : 0);
			// for the variable name:
			datasize += 8 + strlen(temp -> name) + ((strlen(temp -> name) % 8) ? 8 - strlen(temp -> name) % 8 : 0); 
				
			bigbuffer = (int *) calloc(1, datasize);
			offset = 0;
	
			// Always use the miMATRIX header:
			bigbuffer[offset++] = miMATRIX;
			
			// Add up the entire size of the matrix:
			// first the data subelement:
			bigbuffer[offset++] = datasize - 8;
			
			// The array_flags header:
			bigbuffer[offset++] = miUINT32;
			bigbuffer[offset++] = 8;
			// Array_flags values:
			bigbuffer[offset++] = 6 + 8*256*(temp -> type == MATCOMPLEX);
			offset++;
			
			//The dimensions of the matrix:
			bigbuffer[offset++] = miINT32;
			bigbuffer[offset++] = temp -> num_dim * 4;

			// Write the dimensions:
			for(i = 0; i < temp -> num_dim; i++) 
				bigbuffer[i + offset] = temp -> dimensions[i];

			offset += temp -> num_dim;
			if(temp -> num_dim % 2)
				offset ++;

			// The name of the matrix variable:
			bigbuffer[offset++] = miINT8;
			bigbuffer[offset++] = strlen(temp -> name);

			memcpy(bigbuffer + offset, temp -> name, strlen(temp -> name));
			offset += strlen(temp -> name) / 4;
			if(strlen(temp -> name) % 8)
				offset += (strlen(temp -> name) % 8 >= 4) ? 1 : 2;

			// Here we finally write the array data:
			
			if(temp -> type != MATCOMPLEX) {
				bigbuffer[offset++] = miDOUBLE;
				bigbuffer[offset++] = matio_array_size(temp) * 8;
						
				memcpy(bigbuffer + offset, temp -> real, matio_array_size(temp) * 8);
				offset += matio_array_size(temp)*2; // must be multiple of 64 bit since data is double.
			} else {
				bigbuffer[offset++] = miDOUBLE;
				bigbuffer[offset++] = matio_array_size(temp) * 8;
			
				for(i = 0; i < matio_array_size(temp); i ++)
					memcpy(bigbuffer + offset + i * 2, temp -> comp + i, 8);
				
				offset += matio_array_size(temp) * 2; // must be multiple of 64 bit since data is double.
				
				bigbuffer[offset++] = miDOUBLE;
				bigbuffer[offset++] = matio_array_size(temp) * 8;
				
				for(i = 0; i < matio_array_size(temp); i ++)
					memcpy(bigbuffer + offset + i * 2, ((double *)(temp -> comp + i)) + 1, 8);
				
				offset += matio_array_size(temp) * 2;
			}				
#ifdef __ENABLE_WRITE_COMPRESSION	
			compress_size = (unsigned long) (datasize * 1.01 + 12);
			compress_buffer = (unsigned char *) calloc(compress_size, 1);
			status = compress(compress_buffer, &compress_size, (unsigned char *) bigbuffer, datasize);
	
			if(status != Z_OK) // If compression didn't work, just write the uncompressed data.
				fwrite(bigbuffer, datasize, 1, outfile);
			else {
				// First we create the tag indicating that the data is compressed and giving the compressed
				// data size.
				((int *) buffer)[0] = miCOMPRESSED;
				((int *) buffer)[1] = compress_size;
				
				// Then we write this tag, 
				fwrite(buffer, 4, 2, outfile);
				
				// and finally write the compressed buffer.
				fwrite(compress_buffer, compress_size, 1, outfile);
			}
			// Rather inconsistently, compressed arrays within .mat files are not required to be aligned on 64-bit boundaries.
			// We don't have to worry about alignment here.

			free(compress_buffer);
#else
			fwrite(bigbuffer, datasize, 1, outfile);
#endif
			free(bigbuffer);
		}
		temp = temp -> next;
	}
	fclose(outfile);
	return 0;
}

MATdata * matio_complex_MATdata(MATdata *base, int m, int n, char *varname, complex16 *x)
{
	MATdata *ret;
	
	ret = (MATdata *) malloc(sizeof(MATdata));
	ret -> next = base;
	ret -> dimensions = (int *) malloc(2 * sizeof(int));
	ret -> dimensions[0] = m;
	ret -> dimensions[1] = n;
	ret -> num_dim = 2;
	ret -> name = (char *) malloc((strlen(varname) + 1) * sizeof(char));
	strcpy(ret -> name, varname);
	ret -> type = MATCOMPLEX;
	ret -> comp = x;
	ret -> real = NULL;

	return ret;
}

MATdata * matio_real_MATdata(MATdata *base, int m, int n, char *varname, double *x)
{
	MATdata *ret;

	ret = (MATdata *) malloc(sizeof(MATdata));
	ret -> next = base;
	ret -> dimensions = (int *) malloc(2 * sizeof(int));
	ret -> dimensions[0] = m;
	ret -> dimensions[1] = n;
	ret -> num_dim = 2;
	ret -> name = (char *) malloc((strlen(varname) + 1) * sizeof(char));
	strcpy(ret -> name, varname);
	ret -> type = MATREAL;
	ret -> real = x;
	ret -> comp = NULL;

	return ret;
}