function data = loadtiff(fp, img_first, img_last)

% LOADTIFF is a revised version of the tiffread function (v2.0) by
% Francois Nedelec, EMBL (c)1999-2004.
%
% Revision by Xiaolin Nan, UC Berkeley, to integrate into the WFI Reader
% Package.
%
% Last Modified: 10/21/2008

global imgdata;
global h_mainfig;
global TIF;

consolidateStrips = true;

TIF = [];

img_skip  = 0;
img_read  = 0;
img_first = 1;

% set defaults values :

TIF.SampleFormat     = 1;
TIF.SamplesPerPixel  = 1;
TIF.BOS              = 'l';          %byte order string
TIF.file = fp;

% read header
% read byte order: II = little endian, MM = big endian
byte_order = char(fread(TIF.file, 2, 'uchar'));
if ( strcmp(byte_order', 'II') )
    TIF.BOS = 'l';                                %normal PC format
elseif ( strcmp(byte_order','MM') )
    TIF.BOS = 'b';
else
    showmsg(h_mainfig, 'message', 'This is not a TIFF file (no MM or II).');
    data = 0;
    return;
end

%----- read in a number which identifies file as TIFF format
tiff_id = fread(TIF.file,1,'uint16', TIF.BOS);
if (tiff_id ~= 42)
    showmsg(h_mainfig, 'message', 'This is not a TIFF file (missing 42).');
    data = 0;
    return;
end

%----- read the byte offset for the first image file directory (IFD)
ifd_pos = fread(TIF.file,1,'uint32', TIF.BOS);
imgdata = 0;
data = 0;

while (ifd_pos ~= 0)
    clear IMG;

    % IMG.filename = fullfile( pwd, filename );
    % move in the file to the first IFD
    fseek(TIF.file, ifd_pos, -1);

    %disp(strcat('reading img at pos :',num2str(ifd_pos)));
    %read in the number of IFD entries
    num_entries = fread(TIF.file,1,'uint16', TIF.BOS);

    %disp(strcat('num_entries =', num2str(num_entries)));
    %read and process each IFD entry
    for i = 1:num_entries
        % save the current position in the file
        file_pos  = ftell(TIF.file);

        % read entry tag
        TIF.entry_tag = fread(TIF.file, 1, 'uint16', TIF.BOS);
        entry = readIFDentry;
        %disp(strcat('reading entry <',num2str(TIF.entry_tag),'>'));

        switch TIF.entry_tag
            %       case 254
            %           TIF.NewSubfiletype = entry.val;
            case 256         % image width - number of column
                IMG.width          = entry.val;
                data.width         = entry.val;
            case 257         % image height - number of row
                IMG.height         = entry.val;
                TIF.ImageLength    = entry.val;
                data.height        = entry.val;
            case 258         % BitsPerSample per sample
                TIF.BitsPerSample  = entry.val;
                TIF.BytesPerSample = TIF.BitsPerSample / 8;
                IMG.bits           = TIF.BitsPerSample(1);
                %fprintf(1,'BitsPerSample %i %i %i\n', entry.val);
            case 259         % compression
                if (entry.val ~= 1)
                    showmsg(h_mainfig, 'message', 'Compression format not supported.');
                    return;
                end
            case 273         % strip offset
                TIF.StripOffsets   = entry.val;
                TIF.StripNumber    = entry.cnt;
                %fprintf(1,'StripNumber = %i, size(StripOffsets) = %i %i\n', TIF.StripNumber, size(TIF.StripOffsets));
            case 277         % sample_per pixel
                TIF.SamplesPerPixel  = entry.val;
                %fprintf(1,'Color image: sample_per_pixel=%i\n',
                %TIF.SamplesPerPixel);
            case 278         % rows per strip
                TIF.RowsPerStrip   = entry.val;
            case 279         % strip byte counts - number of bytes in each strip after any compressio
                TIF.StripByteCounts= entry.val;
            case 284         %planar configuration describe the order of RGB
                TIF.PlanarConfiguration = entry.val;
                %planar configuration is not fully supported here
                if ( TIF.PlanarConfiguration == 1 )
                    mesg = sprintf('PlanarConfiguration = %i not supported', TIF.PlanarConfiguration);
                    showmsg(h_mainfig, 'message', mesg);
                    data = 0;
                    return;
                end
            case 317        %predictor for compression
                if (entry.val ~= 1)
                    showmsg(h_mainfig, 'message', 'unsuported predictor value');
                    data = 0;
                    return;
                end
            case 339
                TIF.SampleFormat   = entry.val;
                if ( TIF.SampleFormat > 2 )
                    mesg =sprintf('unsuported sample format = %i', TIF.SampleFormat);
                    showmsg(h_mainfig, 'message', mesg);
                    data = 0;
                    return;
                end
            case 33629       %this tag identify the image as a Metamorph stack!
                TIF.MM_stack       = entry.val;
                TIF.MM_stackCnt      = entry.cnt;
        end

        % move to next IFD entry in the file
        fseek(TIF.file, file_pos+12,-1);
    end

    %total number of bytes per image:
    PlaneBytesCnt = IMG.width * IMG.height * TIF.BytesPerSample;

    % prepare a temporary matrix to store normal TIF data
    if ~exist('tempdata', 'var') && ~isfield( TIF, 'MM_stack' )
        % 1.5e8 elements is about the maximum of matrices on a 32-bit XP machine
        max_frames = round(1.2e8 / (IMG.width * IMG.height));
        tempdata = zeros(max_frames, IMG.width, IMG.height, 'uint16');
        %intensity = zeros(max_frames, 1);
    end

    if consolidateStrips
        %Try to consolidate the strips into a single one to speed-up
        %reading:
        BytesCnt = TIF.StripByteCounts(1);

        if BytesCnt < PlaneBytesCnt
            ConsolidateCnt = 1;

            %Count how many Strip are needed to produce a plane
            while TIF.StripOffsets(1) + BytesCnt == TIF.StripOffsets(ConsolidateCnt+1)
                ConsolidateCnt = ConsolidateCnt + 1;
                BytesCnt = BytesCnt + TIF.StripByteCounts(ConsolidateCnt);

                if ( BytesCnt >= PlaneBytesCnt ) 
                    break; 
                end
            end

            %Consolidate the Strips
            if ( BytesCnt <= PlaneBytesCnt(1) ) && ( ConsolidateCnt > 1 )
                %fprintf(1,'Consolidating %i stripes out of %i',
                %ConsolidateCnt, TIF.StripNumber);
                TIF.StripByteCounts = [BytesCnt; TIF.StripByteCounts(ConsolidateCnt+1:TIF.StripNumber ) ];
                TIF.StripOffsets = TIF.StripOffsets( [1 , ConsolidateCnt+1:TIF.StripNumber] );
                TIF.StripNumber  = 1 + TIF.StripNumber - ConsolidateCnt;
            end
        end
    end

    %read the next IFD address:

    ifd_pos = fread(TIF.file, 1, 'uint32', TIF.BOS);
    %if (ifd_pos) disp(['next ifd at', num2str(ifd_pos)]); end

    if isfield( TIF, 'MM_stack' )
        img_last = TIF.MM_stackCnt;

        frames = img_last - img_first + 1;
        imgdata = zeros(frames, width, height, 'uint16');
        intensity = zeros(frames, 1);

        %this loop is to read metamorph stacks:
        % this part will need more testing
        for ii = img_first:img_last
            TIF.StripCnt = 1;

            %read the image
            fileOffset = PlaneBytesCnt * ( ii - 1 );
            %fileOffset = 0;
            %fileOffset = ftell(TIF.file) - TIF.StripOffsets(1);

            if ( TIF.SamplesPerPixel == 1 )
                imgdata(ii, :, :)  = read_plane(fileOffset, IMG.width, IMG.height, 1);
            else
                % IMG.red   = read_plane(fileOffset, IMG.width, IMG.height, 1);
                % IMG.green = read_plane(fileOffset, IMG.width, IMG.height, 2);
                % IMG.blue  = read_plane(fileOffset, IMG.width, IMG.height, 3);
                showmsg('WFI Reader does not handle true color TIF images.');
                data = 0;
                return;
            end

            % print a text timer on the main window, or update the waitbar
            % fprintf(1,'img_read %i img_skip %i\n', img_read, img_skip);
            %[ IMG.info, IMG.MM_stack, IMG.MM_wavelength, IMG.MM_private2 ]
            %= extractMetamorphData(ii);
            img_read = img_read + 1;
            pause(0.001);
            %stack( img_read ) = IMG;
        end
        break;
    else
        %this part to read a normal TIFF stack:
        if ( img_skip + 1 >= img_first )
            TIF.StripCnt = 1;

            %read the image
            if ( TIF.SamplesPerPixel == 1 )
                %fileOffset = PlaneBytesCnt * img_read;
                plane  = read_plane(0, IMG.width, IMG.height, 1);

                if isempty(plane)       % end of file reached.
                    mesg = sprintf('WFI Reader: loading frame %4.0d', img_read);
                    showmsg(h_mainfig, 'message', mesg);
                    break;
                end

                img_read = img_read + 1;
                tempdata(img_read, 1:data.width, 1:data.height) = plane;
                
                if mod(img_read, 10) == 0
                    mesg = sprintf('WFI Reader: loading frame %4.0d', img_read);
                    showmsg(h_mainfig, 'message', mesg);
                end
                %intensity(img_read) = mean(mean(tempdata(img_read, 1:data.width, 1:data.height)));
                pause(0.001);
            else
                %IMG.red   = read_plane(0, IMG.width, IMG.height, 1);
                %IMG.green = read_plane(0, IMG.width, IMG.height, 2);
                %IMG.blue  = read_plane(0, IMG.width, IMG.height, 3);
                showmsg('WFI Reader does not handle true color TIF images.');
                data = 0;
                return;
            end
        else
            img_skip = img_skip + 1;
        end
    end
end

% clean-up and return
imgdata = tempdata(1:img_read, 1:data.width, 1:data.height);
clear tempdata;
clear TIF;
data.frames = img_read;

% compatability assignments
data.badframes = 0;
data.rbad = 0;

if ~ exist( 'data', 'var')
    data = [];
end

return;


%==========================================================================
%==========================================================================
function plane = read_plane(offset, width, height, planeCnt )
global TIF;

%return an empty array if the sample format has zero bits
if ( TIF.BitsPerSample(planeCnt) == 0 )
    plane=[];
    return;
end

%fprintf(1,'reading plane %i size %i %i\n', planeCnt, width, height);
%string description of the type of integer needed: int8 or int16...
typecode = sprintf('int%i', TIF.BitsPerSample(planeCnt) );

%unsigned int if SampleFormat == 1
if ( TIF.SampleFormat == 1 )
    typecode = [ 'u', typecode ];
end

% Preallocate a matrix to hold the sample data:
plane = eval( [ typecode, '(zeros(width, height));'] );

line = 1;
while ( TIF.StripCnt <= TIF.StripNumber )
    strip = read_strip(offset, width, planeCnt, TIF.StripCnt, typecode );

    if isempty(strip)
       plane = [];
       return;
    end

    TIF.StripCnt = TIF.StripCnt + 1;

    % copy the strip onto the data
    plane(:, line:(line+size(strip,2)-1)) = strip;
    line = line + size(strip,2);
    if ( line > height )
        break;
    end
end

% Extract valid part of data if needed
if ~all(size(plane) == [width height]),
    plane = plane(1:width, 1:height);
    error('Cropping data: more bytes read than needed...');
end

% convert the image into 'uint16'
%plane = uint16(plane);
return;

%=================== sub-functions to read a strip ===================
function strip = read_strip(offset, width, planeCnt, stripCnt, typecode)
global TIF;
%fprintf(1,'reading strip at position %i\n',TIF.StripOffsets(stripCnt) +
%offset);
StripLength = TIF.StripByteCounts(stripCnt) ./ TIF.BytesPerSample(planeCnt);
%fprintf(1, 'reading strip %i\n', stripCnt);
fseek(TIF.file, TIF.StripOffsets(stripCnt) + offset, 'bof');
bytes = fread( TIF.file, StripLength, typecode, TIF.BOS );

if ( length(bytes) ~= StripLength )
    strip = [];
    return;
    %error('End of file reached unexpectedly.');
end

bytes = uint16(bytes);
strip = reshape(bytes, width, StripLength / width);
return;

%===================sub-functions that reads an IFD entry:===================
function [nbbytes, typechar] = matlab_type(tiff_typecode)
switch (tiff_typecode)
    case 1
        nbbytes=1;
        typechar='uint8';
    case 2
        nbbytes=1;
        typechar='uchar';
    case 3
        nbbytes=2;
        typechar='uint16';
    case 4
        nbbytes=4;
        typechar='uint32';
    case 5
        nbbytes=8;
        typechar='uint32';
    otherwise
        error('tiff type not supported')
end
return;

%===================sub-functions that reads an IFD entry:===================
function  entry = readIFDentry()
global TIF;

entry.typecode = fread(TIF.file, 1, 'uint16', TIF.BOS);
entry.cnt      = fread(TIF.file, 1, 'uint32', TIF.BOS);
%disp(['typecode =', num2str(entry.typecode),', cnt = ',num2str(entry.cnt)]);

[ entry.nbbytes, entry.typechar ] = matlab_type(entry.typecode);

if entry.nbbytes * entry.cnt > 4
    %next field contains an offset:
    offset = fread(TIF.file, 1, 'uint32', TIF.BOS);

    %disp(strcat('offset = ', num2str(offset)));
    fseek(TIF.file, offset, -1);
end

if TIF.entry_tag == 33629   %special metamorph 'rationals'
    entry.val = fread(TIF.file, 6*entry.cnt, entry.typechar, TIF.BOS);
else
    if entry.typecode == 5
        entry.val = fread(TIF.file, 2*entry.cnt, entry.typechar, TIF.BOS);
    else
        entry.val = fread(TIF.file, entry.cnt, entry.typechar, TIF.BOS);
    end
end

if ( entry.typecode == 2 )
    entry.val = char(entry.val');
end
return;

