/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                  M   M   AAA   TTTTT  L       AAA   BBBB                    %

%                  MM MM  A   A    T    L      A   A  B   B                   %

%                  M M M  AAAAA    T    L      AAAAA  BBBB                    %

%                  M   M  A   A    T    L      A   A  B   B                   %

%                  M   M  A   A    T    LLLLL  A   A  BBBB                    %

%                                                                             %

%                                                                             %

%                        Read MATLAB Image Format                             %

%                                                                             %

%                              Software Design                                %

%                              Jaroslav Fojtik                                %

%                                2001-2008                                    %

%                                                                             %

%                                                                             %

%  Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization         %

%  dedicated to making software imaging solutions freely available.           %

%                                                                             %

%  You may not use this file except in compliance with the License.  You may  %

%  obtain a copy of the License at                                            %

%                                                                             %

%    https://imagemagick.org/license/                                         %

%                                                                             %

%  Unless required by applicable law or agreed to in writing, software        %

%  distributed under the License is distributed on an "AS IS" BASIS,          %

%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %

%  See the License for the specific language governing permissions and        %

%  limitations under the License.                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

*/

/*

  Include declarations.

*/

#include "MagickCore/studio.h"

#include "MagickCore/attribute.h"

#include "MagickCore/blob.h"

#include "MagickCore/blob-private.h"

#include "MagickCore/cache.h"

#include "MagickCore/color-private.h"

#include "MagickCore/colormap.h"

#include "MagickCore/colorspace-private.h"

#include "MagickCore/distort.h"

#include "MagickCore/exception.h"

#include "MagickCore/exception-private.h"

#include "MagickCore/image.h"

#include "MagickCore/image-private.h"

#include "MagickCore/list.h"

#include "MagickCore/magick.h"

#include "MagickCore/memory_.h"

#include "MagickCore/monitor.h"

#include "MagickCore/monitor-private.h"

#include "MagickCore/pixel-accessor.h"

#include "MagickCore/quantum.h"

#include "MagickCore/quantum-private.h"

#include "MagickCore/option.h"

#include "MagickCore/pixel.h"

#include "MagickCore/resource_.h"

#include "MagickCore/static.h"

#include "MagickCore/string_.h"

#include "MagickCore/module.h"

#include "MagickCore/timer-private.h"

#include "MagickCore/transform.h"

#include "MagickCore/utility-private.h"

#include "coders/coders-private.h"

#if defined(MAGICKCORE_ZLIB_DELEGATE)

 #include "zlib.h"

#endif

/*

  Forward declaration.

*/

static MagickBooleanType

  WriteMATImage(const ImageInfo *,Image *,ExceptionInfo *);

/* Auto coloring method, sorry this creates some artefact inside data

MinReal+j*MaxComplex = red  MaxReal+j*MaxComplex = black

MinReal+j*0 = white          MaxReal+j*0 = black

MinReal+j*MinComplex = blue  MaxReal+j*MinComplex = black

*/

typedef struct

{

  char identific[124];

  unsigned short Version;

  char EndianIndicator[2];

  unsigned int DataType;

  unsigned int ObjectSize;

  unsigned int unknown1;

  unsigned int unknown2;

  unsigned short unknown5;

  unsigned char StructureFlag;

  unsigned char StructureClass;

  unsigned int unknown3;

  unsigned int unknown4;

  unsigned int DimFlag;

  unsigned int SizeX;

  unsigned int SizeY;

  unsigned short Flag1;

  unsigned short NameFlag;

}

MATHeader;

static const char

  MonthsTab[12][4] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};

static const char

   DayOfWTab[7][4] = {"Sun","Mon","Tue","Wed","Thu","Fri","Sat"};

static const char

  OsDesc[] =

#if defined(MAGICKCORE_WINDOWS_SUPPORT)

    "PCWIN";

#else

 #ifdef __APPLE__

    "MAC";

 #else

    "LNX86";

 #endif

#endif

typedef enum

  {

    miINT8 = 1,   /* 8 bit signed */

    miUINT8,      /* 8 bit unsigned */

    miINT16,      /* 16 bit signed */

    miUINT16,     /* 16 bit unsigned */

    miINT32,      /* 32 bit signed */

    miUINT32,     /* 32 bit unsigned */

    miSINGLE,     /* IEEE 754 single precision float */

    miRESERVE1,

    miDOUBLE,     /* IEEE 754 double precision float */

    miRESERVE2,

    miRESERVE3,

    miINT64,      /* 64 bit signed */

    miUINT64,     /* 64 bit unsigned */

    miMATRIX,     /* MATLAB array */

    miCOMPRESSED, /* Compressed Data */

    miUTF8,       /* Unicode UTF-8 Encoded Character Data */

    miUTF16,      /* Unicode UTF-16 Encoded Character Data */

    miUTF32       /* Unicode UTF-32 Encoded Character Data */

  } mat5_data_type;

typedef enum

  {

    mxCELL_CLASS=1,  /* cell array */

    mxSTRUCT_CLASS,  /* structure */

    mxOBJECT_CLASS,  /* object */

    mxCHAR_CLASS,    /* character array */

    mxSPARSE_CLASS,  /* sparse array */

    mxDOUBLE_CLASS,  /* double precision array */

    mxSINGLE_CLASS,  /* single precision floating point */

    mxINT8_CLASS,    /* 8 bit signed integer */

    mxUINT8_CLASS,   /* 8 bit unsigned integer */

    mxINT16_CLASS,   /* 16 bit signed integer */

    mxUINT16_CLASS,  /* 16 bit unsigned integer */

    mxINT32_CLASS,   /* 32 bit signed integer */

    mxUINT32_CLASS,  /* 32 bit unsigned integer */

    mxINT64_CLASS,   /* 64 bit signed integer */

    mxUINT64_CLASS,  /* 64 bit unsigned integer */

    mxFUNCTION_CLASS /* Function handle */

  } arrayclasstype;

#define FLAG_COMPLEX 0x8

#define FLAG_GLOBAL  0x4

#define FLAG_LOGICAL 0x2

static const QuantumType z2qtype[4] = {GrayQuantum, BlueQuantum, GreenQuantum, RedQuantum};

static void InsertComplexDoubleRow(Image *image,double *p,int y,double MinVal,

  double MaxVal,ExceptionInfo *exception)

{

  double f;

  int x;

  Quantum *q;

  if (MinVal >= 0)

    MinVal = -1;

  if (MaxVal <= 0)

    MaxVal = 1;

  q=QueueAuthenticPixels(image,0,y,image->columns,1,exception);

  if (q == (Quantum *) NULL)

    return;

  for (x = 0; x < (ssize_t) image->columns; x++)

  {

    if (*p > 0)

      {

        f=(*p/MaxVal)*((double) QuantumRange-(double) GetPixelRed(image,q));

        if (IsNaN(f) != 0)      

          f=0.0;

        if ((f+(double) GetPixelRed(image,q)) >= (double) QuantumRange)

          SetPixelRed(image,QuantumRange,q);

        else

          SetPixelRed(image,GetPixelRed(image,q)+ClampToQuantum(f),q);

        f=(double) GetPixelGreen(image,q)-f/2.0;

        if (IsNaN(f) != 0)      

          f=0.0;

        if (f <= 0.0)

          {

            SetPixelGreen(image,0,q);

            SetPixelBlue(image,0,q);

          }

        else

          {

            SetPixelBlue(image,ClampToQuantum(f),q);

            SetPixelGreen(image,ClampToQuantum(f),q);

          }

      }

    if (*p < 0)

      {

        f=(*p/MinVal)*(Quantum) (QuantumRange-GetPixelBlue(image,q));

        if (IsNaN(f) != 0)      

          f=0.0;

        if ((f+(double) GetPixelBlue(image,q)) >= (double) QuantumRange)

          SetPixelBlue(image,QuantumRange,q);

        else

          SetPixelBlue(image,GetPixelBlue(image,q)+ClampToQuantum(f),q);

        f=GetPixelGreen(image,q)-f/2.0;

        if (f <= 0.0)

          {

            SetPixelRed(image,0,q);

            SetPixelGreen(image,0,q);

          }

        else

          {

            SetPixelRed(image,ClampToQuantum(f),q);

            SetPixelGreen(image,ClampToQuantum(f),q);

          }

      }

    p++;

    q++;

  }

  if (!SyncAuthenticPixels(image,exception))

    return;

  return;

}

static void InsertComplexFloatRow(Image *image,float *p,int y,double MinVal,

  double MaxVal,ExceptionInfo *exception)

{

  double f;

  int x;

  Quantum *q;

  if (MinVal >= 0)

    MinVal = -1;

  if (MaxVal <= 0)

    MaxVal = 1;

  q = QueueAuthenticPixels(image, 0, y, image->columns, 1,exception);

  if (q == (Quantum *) NULL)

    return;

  for (x = 0; x < (ssize_t) image->columns; x++)

  {

    if (*p > 0)

      {

        f=((double) *p/MaxVal)*((double) QuantumRange-(double)

          GetPixelRed(image,q));

        if (IsNaN(f) != 0)      

          f=0.0;

        if ((f+GetPixelRed(image,q)) < QuantumRange)

          SetPixelRed(image,GetPixelRed(image,q)+ClampToQuantum(f),q);

        else

          SetPixelRed(image,QuantumRange,q);

        f/=2.0;

        if (f < GetPixelGreen(image,q))

          {

            SetPixelBlue(image,GetPixelBlue(image,q)-ClampToQuantum(f),q);

            SetPixelGreen(image,GetPixelBlue(image,q),q);

          }

        else

          {

            SetPixelGreen(image,0,q);

            SetPixelBlue(image,0,q);

          }

      }

    if (*p < 0)

      {

        f=(*p/MaxVal)*(Quantum) (QuantumRange-GetPixelBlue(image,q));

        if (IsNaN(f) != 0)      

          f=0.0;

        if ((f+GetPixelBlue(image,q)) < QuantumRange)

          SetPixelBlue(image,GetPixelBlue(image,q)+ClampToQuantum(f),q);

        else

          SetPixelBlue(image,QuantumRange,q);

        f/=2.0;

        if (f < GetPixelGreen(image,q))

          {

            SetPixelRed(image,GetPixelRed(image,q)-ClampToQuantum(f),q);

            SetPixelGreen(image,GetPixelRed(image,q),q);

          }

        else

          {

            SetPixelGreen(image,0,q);

            SetPixelRed(image,0,q);

          }

      }

    p++;

    q++;

  }

  if (!SyncAuthenticPixels(image,exception))

    return;

  return;

}

/************** READERS ******************/

/* This function reads one block of floats*/

static void ReadBlobFloatsLSB(Image * image, size_t len, float *data)

{

  while (len >= 4)

  {

    *data++ = ReadBlobFloat(image);

    len -= sizeof(float);

  }

  if (len > 0)

    (void) SeekBlob(image, (MagickOffsetType) len, SEEK_CUR);

}

static void ReadBlobFloatsMSB(Image * image, size_t len, float *data)

{

  while (len >= 4)

  {

    *data++ = ReadBlobFloat(image);

    len -= sizeof(float);

  }

  if (len > 0)

    (void) SeekBlob(image, (MagickOffsetType) len, SEEK_CUR);

}

/* This function reads one block of doubles*/

static void ReadBlobDoublesLSB(Image * image, size_t len, double *data)

{

  while (len >= 8)

  {

    *data++ = ReadBlobDouble(image);

    len -= sizeof(double);

  }

  if (len > 0)

    (void) SeekBlob(image, (MagickOffsetType) len, SEEK_CUR);

}

static void ReadBlobDoublesMSB(Image * image, size_t len, double *data)

{

  while (len >= 8)

  {

    *data++ = ReadBlobDouble(image);

    len -= sizeof(double);

  }

  if (len > 0)

    (void) SeekBlob(image, (MagickOffsetType) len, SEEK_CUR);

}

/* Calculate minimum and maximum from a given block of data */

static void CalcMinMax(Image *image, int endian_indicator, ssize_t SizeX, ssize_t SizeY, size_t CellType, unsigned ldblk, void *BImgBuff, double *Min, double *Max)

{

MagickOffsetType filepos;

ssize_t i, x;

void (*ReadBlobDoublesXXX)(Image * image, size_t len, double *data);

void (*ReadBlobFloatsXXX)(Image * image, size_t len, float *data);

double *dblrow;

float *fltrow;

  if (endian_indicator == LSBEndian)

  {

    ReadBlobDoublesXXX = ReadBlobDoublesLSB;

    ReadBlobFloatsXXX = ReadBlobFloatsLSB;

  }

  else    /* MI */

  {

    ReadBlobDoublesXXX = ReadBlobDoublesMSB;

    ReadBlobFloatsXXX = ReadBlobFloatsMSB;

  }

  filepos = TellBlob(image);     /* Please note that file seeking occurs only in the case of doubles */

  for (i = 0; i < SizeY; i++)

  {

    if (CellType==miDOUBLE)

    {

      ReadBlobDoublesXXX(image, ldblk, (double *)BImgBuff);

      dblrow = (double *)BImgBuff;

      if (i == 0)

      {

        *Min = *Max = *dblrow;

      }

      for (x = 0; x < SizeX; x++)

      {

        if (*Min > *dblrow)

          *Min = *dblrow;

        if (*Max < *dblrow)

          *Max = *dblrow;

        dblrow++;

      }

    }

    if (CellType==miSINGLE)

    {

      ReadBlobFloatsXXX(image, ldblk, (float *)BImgBuff);

      fltrow = (float *)BImgBuff;

      if (i == 0)

      {

        *Min = *Max = *fltrow;

      }

    for (x = 0; x < (ssize_t) SizeX; x++)

      {

        if (*Min > *fltrow)

          *Min = *fltrow;

        if (*Max < *fltrow)

          *Max = *fltrow;

        fltrow++;

      }

    }

  }

  (void) SeekBlob(image, (MagickOffsetType) filepos, SEEK_SET);

}

static void FixSignedValues(const Image *image,Quantum *q, ssize_t y)

{

  while(y-->0)

  {

     /* Please note that negative values will overflow

        Q=8; QuantumRange=255: <0;127> + 127+1 = <128; 255>

           <-1;-128> + 127+1 = <0; 127> */

    SetPixelRed(image,GetPixelRed(image,q)+QuantumRange/2+1,q);

    SetPixelGreen(image,GetPixelGreen(image,q)+QuantumRange/2+1,q);

    SetPixelBlue(image,GetPixelBlue(image,q)+QuantumRange/2+1,q);

    q++;

  }

}

/** Fix whole row of logical/binary data. It means pack it. */

static void FixLogical(unsigned char *Buff,ssize_t ldblk)

{

unsigned char mask=128;

unsigned char *BuffL = Buff;

unsigned char val = 0;

  if (ldblk == 0)

    return;

  while(ldblk-->0)

  {

    if(*Buff++ != 0)

      val |= mask;

    mask >>= 1;

    if(mask==0)

    {

      *BuffL++ = val;

      val = 0;

      mask = 128;

    }

  }

  *BuffL = val;

}

#if defined(MAGICKCORE_ZLIB_DELEGATE)

static voidpf AcquireZIPMemory(voidpf context,unsigned int items,

  unsigned int size)

{

  (void) context;

  return((voidpf) AcquireQuantumMemory(items,size));

}

static void RelinquishZIPMemory(voidpf context,voidpf memory)

{

  (void) context;

  memory=RelinquishMagickMemory(memory);

}

#endif

#if defined(MAGICKCORE_ZLIB_DELEGATE)

/** This procedure decompreses an image block for a new MATLAB format. */

static Image *decompress_block(Image *orig,unsigned int *Size,

  ImageInfo *clone_info,ExceptionInfo *exception)

{

Image *image2;

void *cache_block, *decompress_block;

z_stream zip_info;

FILE *mat_file;

ssize_t magick_size;

size_t extent;

int file;

MagickBooleanType status;

int zip_status;

ssize_t total_size = 0;

size_t compressed_total = 0;

static const size_t MAX_ABSOLUTE_SIZE = 512UL * 1024UL * 1024UL; /* 512 MB */

static const size_t MAX_EXPANSION_RATIO = 64; /* 64x expansion allowed */

  if (clone_info == NULL)

    return NULL;

  if (clone_info->file)    /* Close file opened from previous transaction. */

    {

      fclose(clone_info->file);

      clone_info->file = NULL;

      (void) remove_utf8(clone_info->filename);

    }

  cache_block=AcquireQuantumMemory((size_t) (*Size < MagickMinBufferExtent) ?

    *Size: MagickMinBufferExtent,sizeof(unsigned char *));

  if (cache_block == NULL)

    return NULL;

  decompress_block=AcquireQuantumMemory((size_t) (4096),

    sizeof(unsigned char *));

  if (decompress_block == NULL)

    {

      RelinquishMagickMemory(cache_block);

      return NULL;

    }

  mat_file=0;

  file=AcquireUniqueFileResource(clone_info->filename);

  if (file != -1)

    mat_file = fdopen(file,"ab+");

  if(!mat_file)

  {

    RelinquishMagickMemory(cache_block);

    RelinquishMagickMemory(decompress_block);

    (void) LogMagickEvent(CoderEvent,GetMagickModule(),"Cannot create file stream for decompressed image");

    return NULL;

  }

  (void) memset(&zip_info,0,sizeof(zip_info));

  zip_info.zalloc=AcquireZIPMemory;

  zip_info.zfree=RelinquishZIPMemory;

  zip_info.opaque=(voidpf) NULL;

  zip_status=inflateInit(&zip_info);

  if (zip_status != Z_OK)

    {

      RelinquishMagickMemory(cache_block);

      RelinquishMagickMemory(decompress_block);

      (void) ThrowMagickException(exception,GetMagickModule(),CorruptImageError,

        "UnableToUncompressImage","`%s'",clone_info->filename);

      (void) fclose(mat_file);

      RelinquishUniqueFileResource(clone_info->filename);

      return(NULL);

    }

  zip_info.avail_in=0;

  zip_info.total_out=0;

  while ((*Size > 0) && (EOFBlob(orig) == 0))

  {

    magick_size=ReadBlob(orig,(*Size < MagickMinBufferExtent) ? *Size :

      MagickMinBufferExtent,(unsigned char *) cache_block);

    if (magick_size == 0)

      break;

    zip_info.next_in=(Bytef *) cache_block;

    zip_info.avail_in=(uInt) magick_size;

    compressed_total+=magick_size;

    while (zip_info.avail_in > 0)

    {

      size_t chunk_size = 0;

      zip_info.avail_out=4096;

      zip_info.next_out=(Bytef *) decompress_block;

      zip_status=inflate(&zip_info,Z_NO_FLUSH);

      if ((zip_status != Z_OK) && (zip_status != Z_STREAM_END))

        break;

      chunk_size=4096-zip_info.avail_out;

      extent=fwrite(decompress_block,1,chunk_size,mat_file);

      if (extent != chunk_size)

        {

          zip_status=Z_DATA_ERROR;

          break;

        }

      total_size+=chunk_size;

      if (compressed_total > 0)

        {

          if ((total_size > (ssize_t) MAX_ABSOLUTE_SIZE) ||

              (total_size > (ssize_t) (compressed_total*MAX_EXPANSION_RATIO)))

            {

              zip_status=Z_DATA_ERROR;

              break;

            }

        }

      if (zip_status == Z_STREAM_END)

        goto DblBreak;

    }

    if ((zip_status != Z_OK) && (zip_status != Z_STREAM_END))

      break;

    *Size-=(unsigned int) magick_size;

  }

DblBreak:

  inflateEnd(&zip_info);

  if (fseek(mat_file,0,SEEK_SET) != 0)

    (void) ThrowMagickException(exception,GetMagickModule(),CorruptImageError,

      "UnableToCreateTemporaryFile","`%s'",clone_info->filename);

  RelinquishMagickMemory(cache_block);

  RelinquishMagickMemory(decompress_block);

  *Size = (unsigned int) total_size;

  clone_info->file=mat_file;

  if( (image2 = AcquireImage(clone_info,exception))==NULL ) goto EraseFile;

  image2->columns=0;

  image2->rows=0;

  status = OpenBlob(clone_info,image2,ReadBinaryBlobMode,exception);

  if (status == MagickFalse)

  {

    DeleteImageFromList(&image2);

EraseFile:

    fclose(clone_info->file);

    clone_info->file = NULL;

    RelinquishUniqueFileResource(clone_info->filename);

    return NULL;

  }

  return(image2);

}

#endif

static Image *ReadMATImageV4(const ImageInfo *image_info,Image *image,

  ExceptionInfo *exception)

{

  typedef struct {

    unsigned char Type[4];

    unsigned int nRows;

    unsigned int nCols;

    unsigned int imagf;

    unsigned int nameLen;

  } MAT4_HDR;

  EndianType

    endian;

  Image

    *rotated_image;

  MagickBooleanType

    status;

  MAT4_HDR

    HDR;

  QuantumInfo

    *quantum_info;

  QuantumFormatType

    format_type;

  ssize_t

    count,

    i,

    ldblk,

    y;

  unsigned char

    *pixels;

  unsigned int

    depth;

  quantum_info=(QuantumInfo *) NULL;

  (void) SeekBlob(image,0,SEEK_SET);

  status=MagickTrue;

  while (EOFBlob(image) == MagickFalse)

  {

    /*

     Object parser loop.

    */

    ldblk=ReadBlobLSBLong(image);

    if(EOFBlob(image)) break;

    if ((ldblk > 9999) || (ldblk < 0))

      break;

    HDR.Type[3]=(unsigned char) (ldblk % 10); ldblk /= 10;  /* T digit */

    HDR.Type[2]=(unsigned char) (ldblk % 10); ldblk /= 10;  /* P digit */

    HDR.Type[1]=(unsigned char) (ldblk % 10); ldblk /= 10;  /* O digit */

    HDR.Type[0]=(unsigned char) (ldblk);        /* M digit */

    if (HDR.Type[3] != 0)

      break;  /* Data format */

    if (HDR.Type[2] != 0)

      break;  /* Always 0 */

    if (HDR.Type[0] == 0)

      {

        HDR.nRows=ReadBlobLSBLong(image);

        HDR.nCols=ReadBlobLSBLong(image);

        HDR.imagf=ReadBlobLSBLong(image);

        HDR.nameLen=ReadBlobLSBLong(image);

        endian=LSBEndian;

      }

    else

      {

        HDR.nRows=ReadBlobMSBLong(image);

        HDR.nCols=ReadBlobMSBLong(image);

        HDR.imagf=ReadBlobMSBLong(image);

        HDR.nameLen=ReadBlobMSBLong(image);

        endian=MSBEndian;

      }

    if ((HDR.imagf != 0) && (HDR.imagf != 1))

      break;

    if (HDR.nameLen > 0xFFFF)

      return(DestroyImageList(image));

    for (i=0; i < (ssize_t) HDR.nameLen; i++)

    {

      int

        byte;

      /*

        Skip matrix name.

      */

      byte=ReadBlobByte(image);

      if (byte == EOF)

        {

          ThrowFileException(exception,CorruptImageError,"UnexpectedEndOfFile",

            image->filename);

          break;

        }

    }

    image->columns=(size_t) HDR.nRows;

    image->rows=(size_t) HDR.nCols;

    if ((image->columns == 0) || (image->rows == 0))

      return(DestroyImageList(image));

    if (image_info->ping != MagickFalse)

      {

        Swap(image->columns,image->rows);

        if(HDR.imagf==1) ldblk *= 2;

        SeekBlob(image, (MagickOffsetType) HDR.nCols*ldblk, SEEK_CUR);

        if ((image->columns == 0) || (image->rows == 0))

          return(image->previous == (Image *) NULL ? DestroyImageList(image)

            : image);

        goto skip_reading_current;

      }

    if ((image->columns > GetBlobSize(image)) ||

        (image->rows > GetBlobSize(image)))

      ThrowReaderException(CorruptImageError,"InsufficientImageDataInFile");

    status=SetImageExtent(image,image->columns,image->rows,exception);

    if (status == MagickFalse)

      return(DestroyImageList(image));

    (void) SetImageBackgroundColor(image,exception);

    (void) SetImageColorspace(image,GRAYColorspace,exception);

    quantum_info=AcquireQuantumInfo(image_info,image);

    if (quantum_info == (QuantumInfo *) NULL)

      return(DestroyImageList(image));

    switch(HDR.Type[1])

    {

      case 0:

        format_type=FloatingPointQuantumFormat;

        depth=64;

        break;

      case 1:

        format_type=FloatingPointQuantumFormat;

        depth=32;

        break;

      case 2:

        format_type=UnsignedQuantumFormat;

        depth=16;

        break;

      case 3:

        format_type=SignedQuantumFormat;

        depth=16;

        break;

      case 4:

        format_type=UnsignedQuantumFormat;

        depth=8;

        break;

      default:

        format_type=UnsignedQuantumFormat;

        depth=8;

        break;

    }

    image->depth=depth;

    if (HDR.Type[0] != 0)

      SetQuantumEndian(image,quantum_info,MSBEndian);

    status=SetQuantumFormat(image,quantum_info,format_type);

    if (status != MagickFalse)

      status=SetQuantumDepth(image,quantum_info,depth);

    if (status != MagickFalse)

      status=SetQuantumEndian(image,quantum_info,endian);

    if (status == MagickFalse)

      {

        quantum_info=DestroyQuantumInfo(quantum_info);

        ThrowReaderException(ResourceLimitError,"MemoryAllocationFailed");

      }

    SetQuantumScale(quantum_info,1.0);

    pixels=(unsigned char *) GetQuantumPixels(quantum_info);

    for (y=0; y < (ssize_t) image->rows; y++)

    {

      Quantum

        *magick_restrict q;

      count=ReadBlob(image,depth/8*image->columns,(char *) pixels);

      if (count == -1)

        break;

      q=QueueAuthenticPixels(image,0,(ssize_t) image->rows-y-1,image->columns,1,

        exception);

      if (q == (Quantum *) NULL)

        break;

      (void) ImportQuantumPixels(image,(CacheView *) NULL,quantum_info,

        GrayQuantum,pixels,exception);

      if ((HDR.Type[1] == 2) || (HDR.Type[1] == 3))

        FixSignedValues(image,q,(int) image->columns);

      if (SyncAuthenticPixels(image,exception) == MagickFalse)

        break;

      if (image->previous == (Image *) NULL)

        {

          status=SetImageProgress(image,LoadImageTag,(MagickOffsetType) y,

            image->rows);

          if (status == MagickFalse)

            break;

        }

    }

    if (HDR.imagf == 1)

      for (y=0; y < (ssize_t) image->rows; y++)

      {

        /*

          Read complex pixels.

        */

        count=ReadBlob(image,depth/8*image->columns,(char *) pixels);

        if (count == -1)

          break;

        if (HDR.Type[1] == 0)

          InsertComplexDoubleRow(image,(double *) pixels,(int) y,0,0,exception);

        else

          InsertComplexFloatRow(image,(float *) pixels,(int) y,0,0,exception);

      }

    quantum_info=DestroyQuantumInfo(quantum_info);

    if (EOFBlob(image) != MagickFalse)

      {

        ThrowFileException(exception,CorruptImageError,"UnexpectedEndOfFile",

          image->filename);

        break;

      }

    rotated_image=RotateImage(image,90.0,exception);

    if (rotated_image != (Image *) NULL)

      {

        rotated_image->page.x=0;

        rotated_image->page.y=0;

        rotated_image->colors = image->colors;

        DestroyBlob(rotated_image);

        rotated_image->blob=ReferenceBlob(image->blob);

        ReplaceImageInList(&image,rotated_image);

        image=rotated_image;

      }

    /*

      Proceed to next image.

    */

    if (image_info->number_scenes != 0)

      if (image->scene >= (image_info->scene+image_info->number_scenes-1))

        break;

    /*

      Allocate next image structure.

    */

skip_reading_current:

    if ((image_info->ping != MagickFalse) && (image_info->number_scenes != 0))

      if (image->scene >= (image_info->scene+image_info->number_scenes-1))

        break;

    AcquireNextImage(image_info,image,exception);

    if (GetNextImageInList(image) == (Image *) NULL)

      {

        status=MagickFalse;

        break;

      }

    image=SyncNextImageInList(image);

    status=SetImageProgress(image,LoadImagesTag,TellBlob(image),

      GetBlobSize(image));

    if (status == MagickFalse)

      break;

  }

  if (CloseBlob(image) == MagickFalse)

    status=MagickFalse;

  if (status == MagickFalse)

    return(DestroyImageList(image));

  return(GetFirstImageInList(image));

}

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

%   R e a d M A T L A B i m a g e                                             %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  ReadMATImage() reads an MAT X image file and returns it.  It

%  allocates the memory necessary for the new Image structure and returns a

%  pointer to the new image.

%

%  The format of the ReadMATImage method is:

%

%      Image *ReadMATImage(const ImageInfo *image_info,ExceptionInfo *exception)

%

%  A description of each parameter follows:

%

%    o image:  Method ReadMATImage returns a pointer to the image after

%      reading. A null image is returned if there is a memory shortage or if

%      the image cannot be read.

%

%    o image_info: Specifies a pointer to a ImageInfo structure.

%

%    o exception: return any errors or warnings in this structure.

%

*/

static Image *ReadMATImage(const ImageInfo *image_info,ExceptionInfo *exception)

{

  Image *image, *image2=NULL,

   *rotated_image;

  Quantum *q;

  unsigned int status;

  MATHeader MATLAB_HDR;

  size_t size;

  size_t CellType;

  QuantumInfo *quantum_info;

  ImageInfo *clone_info;

  ssize_t i;

  ssize_t ldblk;

  unsigned char *BImgBuff = NULL;

  double MinVal, MaxVal;

  unsigned z, z2;

  unsigned Frames;

  MagickBooleanType logging;

  int sample_size;

  MagickOffsetType filepos=0x80;

  unsigned int (*ReadBlobXXXLong)(Image *image);

  unsigned short (*ReadBlobXXXShort)(Image *image);

  void (*ReadBlobDoublesXXX)(Image * image, size_t len, double *data);

  void (*ReadBlobFloatsXXX)(Image * image, size_t len, float *data);

  assert(image_info != (const ImageInfo *) NULL);

  assert(image_info->signature == MagickCoreSignature);

  assert(exception != (ExceptionInfo *) NULL);

  assert(exception->signature == MagickCoreSignature);

  logging = LogMagickEvent(CoderEvent,GetMagickModule(),"enter");

  /*

     Open image file.

   */

  image = AcquireImage(image_info,exception);

  image2 = (Image *) NULL;

  status = OpenBlob(image_info, image, ReadBinaryBlobMode, exception);

  if (status == MagickFalse)

    {

      image=DestroyImageList(image);

      return((Image *) NULL);

    }

  /*

     Read MATLAB image.

   */

  quantum_info=(QuantumInfo *) NULL;

  clone_info=(ImageInfo *) NULL;

  if (ReadBlob(image,124,(unsigned char *) &MATLAB_HDR.identific) != 124)

    ThrowReaderException(CorruptImageError,"ImproperImageHeader");

  if (strncmp(MATLAB_HDR.identific,"MATLAB",6) != 0)

    {

      image=ReadMATImageV4(image_info,image,exception);

      if (image == NULL)

        {

          if ((image != image2) && (image2 != (Image *) NULL))

            image2=DestroyImage(image2);

          if (clone_info != (ImageInfo *) NULL)

            clone_info=DestroyImageInfo(clone_info);

          return((Image *) NULL);

        }

      goto END_OF_READING;

    }

  MATLAB_HDR.Version = ReadBlobLSBShort(image);

  if(ReadBlob(image,2,(unsigned char *) &MATLAB_HDR.EndianIndicator) != 2)

    ThrowReaderException(CorruptImageError,"ImproperImageHeader");

  if (logging)

    (void) LogMagickEvent(CoderEvent,GetMagickModule(),"  Endian %c%c",

      MATLAB_HDR.EndianIndicator[0],MATLAB_HDR.EndianIndicator[1]);

  if (!strncmp(MATLAB_HDR.EndianIndicator, "IM", 2))

  {

    ReadBlobXXXLong = ReadBlobLSBLong;

    ReadBlobXXXShort = ReadBlobLSBShort;

    ReadBlobDoublesXXX = ReadBlobDoublesLSB;

    ReadBlobFloatsXXX = ReadBlobFloatsLSB;

    image->endian = LSBEndian;

  }

  else if (!strncmp(MATLAB_HDR.EndianIndicator, "MI", 2))

  {

    ReadBlobXXXLong = ReadBlobMSBLong;

    ReadBlobXXXShort = ReadBlobMSBShort;

    ReadBlobDoublesXXX = ReadBlobDoublesMSB;

    ReadBlobFloatsXXX = ReadBlobFloatsMSB;

    image->endian = MSBEndian;

  }

  else

    {

MATLAB_KO:

      if ((image != image2) && (image2 != (Image *) NULL))

        image2=DestroyImage(image2);

      if (clone_info != (ImageInfo *) NULL)

        clone_info=DestroyImageInfo(clone_info);

      ThrowReaderException(CorruptImageError,"ImproperImageHeader");

    }

  filepos = TellBlob(image);

  while(filepos < (MagickOffsetType) GetBlobSize(image) && !EOFBlob(image)) /* object parser loop */

  {

    Frames = 1;

    if(filepos > (MagickOffsetType) GetBlobSize(image) || filepos < 0)

      break;

    if(SeekBlob(image,(MagickOffsetType) filepos,SEEK_SET) != filepos) break;

    /* printf("pos=%X\n",TellBlob(image)); */

    MATLAB_HDR.DataType = ReadBlobXXXLong(image);

    if(EOFBlob(image)) break;

    MATLAB_HDR.ObjectSize = ReadBlobXXXLong(image);

    if(EOFBlob(image)) break;

    if((MagickSizeType) (MATLAB_HDR.ObjectSize+filepos) >= GetBlobSize(image))