
// READS OLS FILES FROM THE Olympus confocal microscope
//  should also open most of the TIFF images
//
// Short usage for both OLS and TIFF files
//
//  olsfile="myconfocalimage.ols";
//  //first retrieve the id of the different channels stored in the OLS file
//  ifds=tiff_get_ifds(olsfile);
//  //then we can retrieve the channel we want (here the 1st one)
//  [img,w,h,bps]=tiff_get_image_prop(olsfile,ifds(1));
//  //img is either hypermat w*h*3 for RGB image or mat w*h
//  //bps is bit per sample
//
//  For OLS only:
//
//  meta_data=find_ols_meta_data(olsfile) returns a text with data about the acquisition
//  [xs,ys,zc1,zc2]=get_scale_from_olsfile(olsfile) image dimensions (vectors) in microns

function is_tiff=tiff_check_tiff_file(tiff_file)
    //Check input file really is a tiff file using magic number
    //returns true if file is tiff, false otherwise
    fd=mopen(tiff_file,'rb');
    mseek(2,fd);
    //Tiff file should have a word == 42 at that position
    if 42==mget(1,'s') then
        is_tiff=%t;
    else
        is_tiff=%f;
    end
    mclose(fd);
//    mclose('all');
//    disp('tiff_check_tiff_file');
//    dispfiles();
endfunction

function ifds=tiff_get_ifds(tiff_file)
    //For a valid tiff_file, returns the index of all the IFDs.
    //We assume there is at least one IFD.
    ifds=[];
    fd=mopen(tiff_file,'rb');
    
    ind=4;//where index to first IFD is stored
    mseek(ind,fd);
    IFD=mget(1,'ui');
    ifds=[ifds,IFD];//we assume there is at least one IFD in tiff file
    while IFD~=0
        //number of tags in current IFD
        mseek(IFD,fd);
        n_tags=mget(1,'s');
        //index of next IFD
        mseek(IFD+2+12*n_tags,fd);
        IFD=mget(1,'ui');
        if IFD~=0 then
            ifds=[ifds,IFD];
        end
    end
    mclose(fd);
//    mclose('all');
//    disp('tiff_get_ifds');
//    dispfiles();
endfunction

function [tags_index, tags_id]=tiff_get_ifd_tag(tiff_file,ifd)
    //For a valid tiff_file, and an IFD index, list tags
    //retuns the index of each tag and the id
    tags_index=[];
    tags_id=[];
    fd=mopen(tiff_file,'rb');
        mseek(ifd,fd);
        n_tags=mget(1,'s');
        for i=0:n_tags-1
            tags_index=[tags_index,ifd+2+12*i];
            mseek(ifd+2+12*i,fd);
            tags_id=[tags_id,mget(1,'s')];
        end
    mclose(fd);
//    mclose('all');
//    disp('tiff_get_ifd_tag');
//    dispfiles();
endfunction




function meta_data=find_ols_meta_data(olsfile)
    //tag id code for image description
    _tiff_image_desc=270;
    meta_data=[];
    is_tiff=tiff_check_tiff_file(olsfile);
    if ~is_tiff then
        error(olsfile+" does not seem to be a valid OLS or TIFF file. Giving up");
        return;
    end
    ifds=tiff_get_ifds(olsfile)
//    fd=mopen(olsfile,'r');
    for ifd=ifds
        [tags_index, tags_id]=tiff_get_ifd_tag(olsfile,ifd);
        rg=find(tags_id==_tiff_image_desc);
        fd=mopen(olsfile,'rb');
        if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);
            
            mseek(tags_index(rg)+8,fd);
            offset=mget(1,'ui',fd);
            
            mseek(offset,fd);
            charcode=mget(n_val,'c',fd);
            str=ascii2string(charcode);
            meta_data=[meta_data;tokens(str, ascii(10))];
        end
        mclose(fd);
    end
//    mclose(fd);
//    mclose('all');
//    disp('find_ols_meta_data');
//    dispfiles();
endfunction





function [img,w,h,bps]=OLD_tiff_get_image_prop(tiff_file,ifd)
    //For a valid tiff_file, and an IFD index
    //
    //return image properties
    //
    _tiff_image_desc=270;//tag id code for image description
    _tiff_image_width=256;//tag id code for image width
    _tiff_image_length=257;//tag id code for image Length
    _tiff_image_bps=258;//tag id code for image bps
    _tiff_image_comp=259;//tag id code for image compression
    _tiff_image_phot_int=262;//tag id code for photometric interpretation
    _tiff_image_stripoffset=273;//tag id code for strip offset
    _tiff_image_orientation=274;//tag id code for image orientation
    _tiff_image_spp=277;//tag id code for samples per pixel
    _tiff_image_rps=278;//tag id code for row per strip
    _tiff_image_sbc=279;//tag id code for strip byte count
    _tiff_xres=282;//tag id code for x resolution
    _tiff_xres=283;//tag id code for y resolution
    _tiff_planarconf=284;//tag id code for planar configuration
    _tiff_resunit=296;//tag id code for resolution unit
    _tiff_date=306;//tag id code for date/time
    _tiff_cmap=320;//tag id code for colormap

    w=[];
    h=[];
    bps=[];
    spp=[];
    comp_=[];
    [tags_index, tags_id]=tiff_get_ifd_tag(tiff_file,ifd);
    
    // I should use datatype to set mget to retrieve the right format
    //  1   -> byte=word=short?
    //  2   -> ascii=char
    //  4   -> long=unsigned integer
    //  5   -> rational ie 2 long that defines ratio long1/long2
    dtype=['s'; 'c'; 'never_here';'ui'; 'man_what_the_hell'];
    
    
    fd=mopen(tiff_file,'rb');    
    rg=find(tags_id==_tiff_image_comp);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        comp_=mget(1,'s',fd);
     end
     
     if comp_~=1 then
         warning("comp_="+string(comp_));
//         error("Cannot handle compressed Tiff");
//         return;
     end
    
    rg=find(tags_id==_tiff_image_phot_int);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        pi_=mget(1,'s',fd);
     end
     
     if pi_==1 then
         disp("grayscale image");
     else
         disp("pi_="+string(pi_));
     end
     
        
    rg=find(tags_id==_tiff_image_bps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        if n_val==1 then
            bps=mget(n_val,'s',fd);
        else
            offset=mget(1,'ui',fd);
            mseek(offset,fd);
            bps=mget(n_val,'s',fd);
        end
     end
     //let's be lazy for sample per pixel
     spp=length(bps);
//    pause 
    rg=find(tags_id==_tiff_image_length);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        h=mget(1,'s',fd);
     end
 
    rg=find(tags_id==_tiff_image_width);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        w=mget(1,'s',fd);
     end
    
//    pause
    rg=find(tags_id==_tiff_image_stripoffset);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        strips_per_image_=n_val;
        mseek(tags_index(rg)+8,fd);
        strip_offset_=mget(1,dtype(data_type),fd);
        disp("strips_per_image_="+string(strips_per_image_));
        disp("strip_offset_="+string(strip_offset_));
     end
     
     //handle row per strip
    rg=find(tags_id==_tiff_image_rps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        pause
        row_per_strip_=mget(1,dtype(data_type),fd);
        disp(" row_per_strip_="+string( row_per_strip_));
        disp("n_val="+string(n_val));
     end
     
     
     rg=find(tags_id==_tiff_image_sbc);
     if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        strip_byte_counts_=mget(1,dtype(data_type),fd);
        disp("strip_byte_counts_="+string(strip_byte_counts_));
        disp("n_val="+string(n_val));
     end
     
     //we go where image starts
     mseek(strip_offset_,fd);
     bytesleft=strip_byte_counts_;
//     pause
     if length(bps)==3 & bps(1)==8 then
        //we have 8 bit rgb
        r=uint8(zeros(h*w,1));g=r;b=r;

        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        
        //both should be the same
        disp(strip_byte_counts_);
        disp(size(tmp,2));
        //hack for windows
        if size(tmp,2)~=strip_byte_counts_  then
            fd1=mopen(tiff_file);
            mseek(strip_offset_,fd1);
            tmp1=mgeti(strip_byte_counts_,'ucl',fd1);
            mclose(fd1);
            disp("Windows sucks")
            disp(strip_byte_counts_);
            disp(size(tmp1,2));
            tmp=tmp1;
//            pause
        else
            disp('Everything is fine under Linux');
        end

        
        r=tmp(1:3:$);
        g=tmp(2:3:$);
        b=tmp(3:3:$);
        r=matrix(r,w,h)';
        g=matrix(g,w,h)';
        b=matrix(b,w,h)';
        img(:,:,1)=r;
        img(:,:,2)=g;
        img(:,:,3)=b;
//        pause
    elseif size(bps)==1 & bps(1)==16 then
        //we have 16bit grayscale
        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    elseif size(bps)==1 & bps(1)==8 then
        //we have 8bit grayscale
        // this is not OLS
        disp("8bit grayscale")
        disp("this does not come from OLS file, I will try to handle it");
//        pause
        img=[];
//        TODO: implement it!
//  http://www.fileformat.info/format/tiff/egff.htm
// I think that for example here:
//  fname='/usr/local/PHOTO/Techno/Bati_Epi/Riber412/A-76/A76-B4--test-de-dose/MEB-ap-HF/A76-B4-apHF-curved-elements-D_q07.tif';
//  ifds=tiff_get_ifds(fname);
//  [img,w,h,bps]=tiff_get_image_prop(fname,ifds(1));
// the tag containing strip bytes counts contains 640 strips and it points to an array of 640 values that is the number of bytes for each strip
//  see explanation here for tag data that are too long to fit in the tag: http://www.awaresystems.be/imaging/tiff/faq.html
        rg=find(tags_id==_tiff_image_sbc);
            if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);//if n_val>1 we have a multistrip tiff image
            mseek(tags_index(rg)+8,fd);
            strip_byte_counts_=mget(1,dtype(data_type),fd);
            mseek(strip_byte_counts_,fd);//strip_byte_counts_ is just a pointer if n_Val>1 as data does not fit in
            //HOURRA the number of bytes for each strip
            sbc_array_=mget(n_val,dtype(data_type),fd);
        end
        //we need the offset for each strip
        rg=find(tags_id==_tiff_image_stripoffset);
        if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);
            strips_per_image_=n_val;
            mseek(tags_index(rg)+8,fd);
            strip_offset_=mget(1,dtype(data_type),fd);
            mseek(strip_offset_,fd);//strip_offset_ is just a pointer if n_Val>1 as data does not fit in
            soffset_array_=mget(n_val,dtype(data_type),fd);
         end
         
//        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'uc',fd);
//        pause
        img=matrix(tmp,w,h).';
//        
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)'; 
//         
//         //now we can read:
//         tmp=[];
//         for i=1:length(soffset_array_)
//             mseek(soffset_array_(i));
//             tmp=[tmp,mget(sbc_array_(i),'uc',fd)];
//         end
//         img=uint8(matrix(tmp,w,h))';
////        pause
////        mseek(strip_offset_,fd);
////        tmp=mgeti(strip_byte_counts_,'uc',fd);
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    end
     mclose(fd);
endfunction


function [img,w,h,bps]=tiff_get_image_prop(tiff_file,ifd)
    //For a valid tiff_file, and an IFD index
    //
    //return image properties
    //
    _tiff_image_desc=270;//tag id code for image description
    _tiff_image_width=256;//tag id code for image width
    _tiff_image_length=257;//tag id code for image Length
    _tiff_image_bps=258;//tag id code for image bps
    _tiff_image_comp=259;//tag id code for image compression
    _tiff_image_phot_int=262;//tag id code for photometric interpretation
    _tiff_image_stripoffset=273;//tag id code for strip offset
    _tiff_image_orientation=274;//tag id code for image orientation
    _tiff_image_spp=277;//tag id code for samples per pixel
    _tiff_image_rps=278;//tag id code for row per strip
    _tiff_image_sbc=279;//tag id code for strip byte count
    _tiff_xres=282;//tag id code for x resolution
    _tiff_xres=283;//tag id code for y resolution
    _tiff_planarconf=284;//tag id code for planar configuration
    _tiff_resunit=296;//tag id code for resolution unit
    _tiff_date=306;//tag id code for date/time
    _tiff_cmap=320;//tag id code for colormap

    w=[];
    h=[];
    bps=[];
    spp=[];
    comp_=[];
    [tags_index, tags_id]=tiff_get_ifd_tag(tiff_file,ifd);
    
    // I should use datatype to set mget to retrieve the right format
    //  1   -> byte=word=short?
    //  2   -> ascii=char
    //  4   -> long=unsigned integer
    //  5   -> rational ie 2 long that defines ratio long1/long2
    dtype=['uc'; 'c'; 'us';'ui'; 'rational'; 'c'; 'undefine'; 's'; 'i'; 'srational'; 'f'; 'd' ];
    
    
    fd=mopen(tiff_file,'rb');    
    rg=find(tags_id==_tiff_image_comp);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        comp_=mget(1,'s',fd);
     end
     
     if comp_~=1 then // it seems that 32771 is also valid for uncompressed
         warning("comp_="+string(comp_));
         if comp_==5 then // from http://www.fileformat.info/format/tiff/egff.htm#X058-9-TIFF-TB-3
            warning("Looks k=like LZW compression");
            // we could try to implement it: https://www.fileformat.info/format/tiff/corion-lzw.htm
            // https://rosettacode.org/wiki/LZW_compression#Julia
            // each strip is compressed separately: comparing compressed/uncompressed might help.
         end
         
         error("Cannot handle compressed Tiff");
         return;
     end
    
    rg=find(tags_id==_tiff_image_phot_int);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        pi_=mget(1,'s',fd);
     end
     
     if pi_==1 then
         disp("grayscale image");
     else
         disp("pi_="+string(pi_));
     end
     
        
    rg=find(tags_id==_tiff_image_bps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        if n_val==1 then
            bps=mget(n_val,'s',fd);
        else
            offset=mget(1,'ui',fd);
            mseek(offset,fd);
            bps=mget(n_val,'s',fd);
        end
     end
     //let's be lazy for sample per pixel
     spp=length(bps);
//    pause 
    rg=find(tags_id==_tiff_image_length);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        h=mget(1,'s',fd);
     end
 
    rg=find(tags_id==_tiff_image_width);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        w=mget(1,'s',fd);
     end
    
//    pause
    rg=find(tags_id==_tiff_image_stripoffset);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        strips_per_image_=n_val;
        mseek(tags_index(rg)+8,fd);
        strip_offset_=mget(1,dtype(data_type),fd);
        disp("strips_per_image_="+string(strips_per_image_));
        disp("strip_offset_="+string(strip_offset_));
     end
     
     //handle row per strip
    rg=find(tags_id==_tiff_image_rps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
//        pause
        row_per_strip_=mget(1,dtype(data_type),fd);
        disp(" row_per_strip_="+string( row_per_strip_));
        disp("n_val="+string(n_val));
        //here we are: each strip contains row_per_strip_/w  lines (am I transposing something?)
     end
     
     
     rg=find(tags_id==_tiff_image_sbc);
     if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        strip_byte_counts_=mget(1,dtype(data_type),fd);
        disp("strip_byte_counts_="+string(strip_byte_counts_));
        disp("n_val="+string(n_val));
     end
     
     //we go where image starts
     mseek(strip_offset_,fd);
     bytesleft=strip_byte_counts_;
//     pause
     if length(bps)==3 & bps(1)==8 then
        //we have 8 bit rgb
        r=uint8(zeros(h*w,1));g=r;b=r;

        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        
        //both should be the same
        disp(strip_byte_counts_);
        disp(size(tmp,2));
        //hack for windows
        if size(tmp,2)~=strip_byte_counts_  then
            fd1=mopen(tiff_file);
            mseek(strip_offset_,fd1);
            tmp1=mgeti(strip_byte_counts_,'ucl',fd1);
            mclose(fd1);
            disp("Windows sucks")
            disp(strip_byte_counts_);
            disp(size(tmp1,2));
            tmp=tmp1;
//            pause
        else
            disp('Everything is fine under Linux');
        end

        
        r=tmp(1:3:$);
        g=tmp(2:3:$);
        b=tmp(3:3:$);
        r=matrix(r,w,h)';
        g=matrix(g,w,h)';
        b=matrix(b,w,h)';
        img(:,:,1)=r;
        img(:,:,2)=g;
        img(:,:,3)=b;
//        pause
    elseif size(bps)==1 & bps(1)==16 then
        //we have 16bit grayscale
        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    elseif size(bps)==1 & bps(1)==8 then
        //we have 8bit grayscale
        // this is not OLS
        disp("8bit grayscale")
        disp("this does not come from OLS file, I will try to handle it");
//        pause
        img=[];
//        TODO: implement it!
//  http://www.fileformat.info/format/tiff/egff.htm
// I think that for example here:
//  fname='/usr/local/PHOTO/Techno/Bati_Epi/Riber412/A-76/A76-B4--test-de-dose/MEB-ap-HF/A76-B4-apHF-curved-elements-D_q07.tif';
//  ifds=tiff_get_ifds(fname);
//  [img,w,h,bps]=tiff_get_image_prop(fname,ifds(1));
// the tag containing strip bytes counts contains 640 strips and it points to an array of 640 values that is the number of bytes for each strip
//  see explanation here for tag data that are too long to fit in the tag: http://www.awaresystems.be/imaging/tiff/faq.html
        rg=find(tags_id==_tiff_image_sbc);
            if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);//if n_val>1 we have a multistrip tiff image
            mseek(tags_index(rg)+8,fd);
            strip_byte_counts_=mget(1,dtype(data_type),fd);
            mseek(strip_byte_counts_,fd);//strip_byte_counts_ is just a pointer if n_Val>1 as data does not fit in
            //HOURRA the number of bytes for each strip
            sbc_array_=mget(n_val,dtype(data_type),fd);
        end
        //we need the offset for each strip
        rg=find(tags_id==_tiff_image_stripoffset);
        if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);
            strips_per_image_=n_val;
            mseek(tags_index(rg)+8,fd);
            strip_offset_=mget(1,dtype(data_type),fd);
            mseek(strip_offset_,fd);//strip_offset_ is just a pointer if n_Val>1 as data does not fit in
            soffset_array_=mget(n_val,dtype(data_type),fd);
         end
         
//        mseek(strip_offset_,fd);
//      we need to read sbc_array_ byte at each soffset_array_ but we don't read a byte at a time, it depend on the bps for bps=8 uc is fine
        if bps==8 then
            tmp=[];//would be better to allocate the right size to start with
            for i=1:length(soffset_array_)
                curr_offset=soffset_array_(i);
                curr_nbyte_to_read=sbc_array_(i);
                 mseek(curr_offset,fd);
                tmp=[tmp,mgeti(curr_nbyte_to_read,'uc',fd)];
            end
        end
    
//        pause
        img=matrix(tmp,w,h).';
//        
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)'; 
//         
//         //now we can read:
//         tmp=[];
//         for i=1:length(soffset_array_)
//             mseek(soffset_array_(i));
//             tmp=[tmp,mget(sbc_array_(i),'uc',fd)];
//         end
//         img=uint8(matrix(tmp,w,h))';
////        pause
////        mseek(strip_offset_,fd);
////        tmp=mgeti(strip_byte_counts_,'uc',fd);
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    end
     mclose(fd);
endfunction

function [xs,ys,zc1,zc2]=get_scale_from_olsfile(olsfile)
    meta_data=find_ols_meta_data(olsfile);
    execstr(strsubst(meta_data(grep(meta_data,'Y Coordinate')),'Y Coordinate','y0'));
    y0=y0/1000;//y stage position in micron
    execstr(strsubst(meta_data(grep(meta_data,'X Coordinate')),'X Coordinate','x0'));
    x0=x0/1000;//x stage position in micron
    
    ids=grep(meta_data,'[Data 1 Info]')
    ide=grep(meta_data,'[Data 1 Info End]')
    
    data1_info=meta_data(ids:ide);
    ind=grep(data1_info,'XY Convert Value');
    execstr(strsubst(data1_info(ind),'XY Convert Value','xyc1'));
    xyc1=xyc1/1000;//xy pixel size in micron
    
    ind=grep(data1_info,'Z Convert Value');
    execstr(strsubst(data1_info(ind),'Z Convert Value','zc1'));
    zc1=zc1/1000;//z scale unit size in micron
    
    ids=grep(meta_data,'[Data 2 Info]')
    ide=grep(meta_data,'[Data 2 Info End]')
    
    data2_info=meta_data(ids:ide);
    ind=grep(data2_info,'XY Convert Value');
    execstr(strsubst(data2_info(ind),'XY Convert Value','xyc2'));
    xyc2=xyc2/1000;//xy pixel size in micron
    
    ind=grep(data2_info,'Z Convert Value');
    execstr(strsubst(data2_info(ind),'Z Convert Value','zc2'));
    zc2=zc2/1000;//z scale unit size in micron
    
    //dirty dirty: magic constants
    xs=x0+xyc1*([1:1024]-1024/2);
    ys=y0-xyc1*([1:768]-768/2);
endfunction


function [img,w,h,bps]=ReadTiffImage(tiff_file)
    is_tiff=tiff_check_tiff_file(tiff_file);
    if ~is_tiff then
        error("This is not a tiff file.")
    else
        ifds=tiff_get_ifds(tiff_file);
        if isempty(ifds) then
            error("Cannot find a image file directory inside the tiff file")
        else
            //we read the first image we found
            [img,w,h,bps]=tiff_get_image_prop(tiff_file,ifds(1));
        end
        
    end
endfunction

function [strips]=tiff_get_image_strips(tiff_file,ifd)
    //For a valid tiff_file, and an IFD index
    //
    //return image properties
    //
    _tiff_image_desc=270;//tag id code for image description
    _tiff_image_width=256;//tag id code for image width
    _tiff_image_length=257;//tag id code for image Length
    _tiff_image_bps=258;//tag id code for image bps
    _tiff_image_comp=259;//tag id code for image compression
    _tiff_image_phot_int=262;//tag id code for photometric interpretation
    _tiff_image_stripoffset=273;//tag id code for strip offset
    _tiff_image_orientation=274;//tag id code for image orientation
    _tiff_image_spp=277;//tag id code for samples per pixel
    _tiff_image_rps=278;//tag id code for row per strip
    _tiff_image_sbc=279;//tag id code for strip byte count
    _tiff_xres=282;//tag id code for x resolution
    _tiff_xres=283;//tag id code for y resolution
    _tiff_planarconf=284;//tag id code for planar configuration
    _tiff_resunit=296;//tag id code for resolution unit
    _tiff_date=306;//tag id code for date/time
    _tiff_cmap=320;//tag id code for colormap

    w=[];
    h=[];
    bps=[];
    spp=[];
    comp_=[];
    [tags_index, tags_id]=tiff_get_ifd_tag(tiff_file,ifd);
    
    // I should use datatype to set mget to retrieve the right format
    //  1   -> byte=word=short?
    //  2   -> ascii=char
    //  4   -> long=unsigned integer
    //  5   -> rational ie 2 long that defines ratio long1/long2
    dtype=['uc'; 'c'; 'us';'ui'; 'rational'; 'c'; 'undefine'; 's'; 'i'; 'srational'; 'f'; 'd' ];
    
    
    fd=mopen(tiff_file,'rb');    
    rg=find(tags_id==_tiff_image_comp);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        comp_=mget(1,'s',fd);
     end
     
     if comp_~=1 then // it seems that 32771 is also valid for uncompressed
         warning("comp_="+string(comp_));
         if comp_==5 then // from http://www.fileformat.info/format/tiff/egff.htm#X058-9-TIFF-TB-3
            warning("Looks k=like LZW compression");
            // we could try to implement it: https://www.fileformat.info/format/tiff/corion-lzw.htm
            // https://rosettacode.org/wiki/LZW_compression#Julia
            // each strip is compressed separately: comparing compressed/uncompressed might help.
         end
         
//         error("Cannot handle compressed Tiff");
//         return;
     end
    
    rg=find(tags_id==_tiff_image_phot_int);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        pi_=mget(1,'s',fd);
     end
     
     if pi_==1 then
         disp("grayscale image");
     else
         disp("pi_="+string(pi_));
     end
     
        
    rg=find(tags_id==_tiff_image_bps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        if n_val==1 then
            bps=mget(n_val,'s',fd);
        else
            offset=mget(1,'ui',fd);
            mseek(offset,fd);
            bps=mget(n_val,'s',fd);
        end
     end
     //let's be lazy for sample per pixel
     spp=length(bps);
//    pause 
    rg=find(tags_id==_tiff_image_length);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        h=mget(1,'s',fd);
     end
 
    rg=find(tags_id==_tiff_image_width);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        
        mseek(tags_index(rg)+8,fd);
        w=mget(1,'s',fd);
     end
    
//    pause
    rg=find(tags_id==_tiff_image_stripoffset);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        strips_per_image_=n_val;
        mseek(tags_index(rg)+8,fd);
        strip_offset_=mget(1,dtype(data_type),fd);
        disp("strips_per_image_="+string(strips_per_image_));
        disp("strip_offset_="+string(strip_offset_));
     end
     
     //handle row per strip
    rg=find(tags_id==_tiff_image_rps);
    if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
//        pause
        row_per_strip_=mget(1,dtype(data_type),fd);
        disp(" row_per_strip_="+string( row_per_strip_));
        disp("n_val="+string(n_val));
        //here we are: each strip contains row_per_strip_/w  lines (am I transposing something?)
     end
     
     
     rg=find(tags_id==_tiff_image_sbc);
     if rg~=[] then
        mseek(tags_index(rg)+2,fd);
        data_type=mget(1,'s',fd);
        mseek(tags_index(rg)+4,fd);
        n_val=mget(1,'ui',fd);
        mseek(tags_index(rg)+8,fd);
        strip_byte_counts_=mget(1,dtype(data_type),fd);
        disp("strip_byte_counts_="+string(strip_byte_counts_));
        disp("n_val="+string(n_val));
     end
     
     //we go where image starts
     mseek(strip_offset_,fd);
     bytesleft=strip_byte_counts_;
//     pause
     if length(bps)==3 & bps(1)==8 then
        //we have 8 bit rgb
        r=uint8(zeros(h*w,1));g=r;b=r;

        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        
        //both should be the same
        disp(strip_byte_counts_);
        disp(size(tmp,2));
        //hack for windows
        if size(tmp,2)~=strip_byte_counts_  then
            fd1=mopen(tiff_file);
            mseek(strip_offset_,fd1);
            tmp1=mgeti(strip_byte_counts_,'ucl',fd1);
            mclose(fd1);
            disp("Windows sucks")
            disp(strip_byte_counts_);
            disp(size(tmp1,2));
            tmp=tmp1;
//            pause
        else
            disp('Everything is fine under Linux');
        end

        
        r=tmp(1:3:$);
        g=tmp(2:3:$);
        b=tmp(3:3:$);
        r=matrix(r,w,h)';
        g=matrix(g,w,h)';
        b=matrix(b,w,h)';
        img(:,:,1)=r;
        img(:,:,2)=g;
        img(:,:,3)=b;
//        pause
    elseif size(bps)==1 & bps(1)==16 then
        //we have 16bit grayscale
        mseek(strip_offset_,fd);
        tmp=mgeti(strip_byte_counts_,'ucl',fd);
        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    elseif size(bps)==1 & bps(1)==8 then
        //we have 8bit grayscale
        // this is not OLS
        disp("8bit grayscale")
        disp("this does not come from OLS file, I will try to handle it");
//        pause
        img=[];
//        TODO: implement it!
//  http://www.fileformat.info/format/tiff/egff.htm
// I think that for example here:
//  fname='/usr/local/PHOTO/Techno/Bati_Epi/Riber412/A-76/A76-B4--test-de-dose/MEB-ap-HF/A76-B4-apHF-curved-elements-D_q07.tif';
//  ifds=tiff_get_ifds(fname);
//  [img,w,h,bps]=tiff_get_image_prop(fname,ifds(1));
// the tag containing strip bytes counts contains 640 strips and it points to an array of 640 values that is the number of bytes for each strip
//  see explanation here for tag data that are too long to fit in the tag: http://www.awaresystems.be/imaging/tiff/faq.html
        rg=find(tags_id==_tiff_image_sbc);
            if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);//if n_val>1 we have a multistrip tiff image
            mseek(tags_index(rg)+8,fd);
            strip_byte_counts_=mget(1,dtype(data_type),fd);
            mseek(strip_byte_counts_,fd);//strip_byte_counts_ is just a pointer if n_Val>1 as data does not fit in
            //HOURRA the number of bytes for each strip
            sbc_array_=mget(n_val,dtype(data_type),fd);
        end
        //we need the offset for each strip
        rg=find(tags_id==_tiff_image_stripoffset);
        if rg~=[] then
            mseek(tags_index(rg)+2,fd);
            data_type=mget(1,'s',fd);
            mseek(tags_index(rg)+4,fd);
            n_val=mget(1,'ui',fd);
            strips_per_image_=n_val;
            mseek(tags_index(rg)+8,fd);
            strip_offset_=mget(1,dtype(data_type),fd);
            mseek(strip_offset_,fd);//strip_offset_ is just a pointer if n_Val>1 as data does not fit in
            soffset_array_=mget(n_val,dtype(data_type),fd);
         end
         
//        mseek(strip_offset_,fd);
//      we need to read sbc_array_ byte at each soffset_array_ but we don't read a byte at a time, it depend on the bps for bps=8 uc is fine
        if bps==8 then
//            pause
            if n_val~=1 then // we have multiple stripes
                strips=cell(1,length(soffset_array_));
                for i=1:length(soffset_array_)
                    curr_offset=soffset_array_(i);
                    curr_nbyte_to_read=sbc_array_(i);
                     mseek(curr_offset,fd);
                    strips{i}=mgeti(curr_nbyte_to_read,'uc',fd);
                end
            else // one single stripe
                strips=cell(1,1);
                mseek(strip_offset_,fd);
                strips{1}=mgeti(strip_byte_counts_,'uc',fd);
            end
        end
    
//        pause
//        img=matrix(tmp,w,h).';
//        
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)'; 
//         
//         //now we can read:
//         tmp=[];
//         for i=1:length(soffset_array_)
//             mseek(soffset_array_(i));
//             tmp=[tmp,mget(sbc_array_(i),'uc',fd)];
//         end
//         img=uint8(matrix(tmp,w,h))';
////        pause
////        mseek(strip_offset_,fd);
////        tmp=mgeti(strip_byte_counts_,'uc',fd);
////        img=matrix(double(tmp(1:2:$))+2^8*double(tmp(2:2:$)),w,h)';
    end
     mclose(fd);
endfunction

//////
////////olsfiles=listfiles('/host/Dropbox/ap-enduction-lext/2107-28-ap-enduction004_*.ols');
////////olsfiles=listfiles('~/Dropbox/ap-enduction-lext/2107-28-ap-enduction004_*.ols');
////////
//////
////////disp("hep");
////////mclose('all');
//////
////////olsfiles=listfiles('Y:/Techno/Plaques LAAS/2107/2107-28/ap-enduction-lext/2107-28-ap-enduction004_*.ols');
//////olsfiles=listfiles('/usr/local/diodesLASER/Techno/Plaques LAAS/2107/2107-28/ap-enduction-lext/2107-28-ap-enduction004_*.ols');
////////h=scf(0);
////////h.color_map=hotcolormap(16);
//////
//////nf=length(length(olsfiles));
//////
//////allx=[];
//////ally=[];
//////allimg=[];
//////for i=1:nf do
////////    i=13
////////    clf(0);
//////    olsfile=olsfiles(i);
////////    ifds=tiff_get_ifds(olsfile);
////////    [img,w,h,bps]=tiff_get_image_prop(olsfile,ifds(1));
////////    imshow(img)
//////    [xs,ys,zc1,zc2]=get_scale_from_olsfile(olsfile);
//////    allx=[allx,xs];
//////    ally=[ally,ys];
////////    Sgrayplot(xs,ys,double(img(:,:,1))'*zc1);
////////    pause
////////    plot(xs,ones(xs)*min(ys),'.-');
////////    plot(xs,ones(xs)*max(ys),'.-');
////////    plot(ones(ys)*min(xs),ys,'.-');
////////    plot(ones(ys)*max(xs),ys,'.-');
//////end
//////
//////
//////xx=gsort(unique(allx),'g','i');
//////yy=gsort(unique(ally),'g','i');
//////
//////bigimg=uint8(ones(yy'*xx));//only one canal only for the moment
//////bigimgr=uint8(ones(yy'*xx));//only one canal only for the moment
//////bigimgg=uint8(ones(yy'*xx));//only one canal only for the moment
//////bigimgb=uint8(ones(yy'*xx));//only one canal only for the moment
//////
//////for i=1:nf do
////////    pause
//////    olsfile=olsfiles(i);
//////    ifds=tiff_get_ifds(olsfile);
//////    [img,w,h,bps]=tiff_get_image_prop(olsfile,ifds(1));
//////    [xs,ys,zc1,zc2]=get_scale_from_olsfile(olsfile);
//////    rgx=find(xx>=min(xs) & xx<=max(xs));
//////    rgy=find(yy>=min(ys) & yy<=max(ys));
//////    yyM=max(yy);
//////    if size(img,1)==size(rgy,2) & size(img,2)==size(rgx,2) then
//////        //perfect match pixel by pixel
//////        bigimg(rgy($:-1:1),rgx)=img(:,:,1);
//////        bigimgr(rgy($:-1:1),rgx)=img(:,:,1);
//////        bigimgg(rgy($:-1:1),rgx)=img(:,:,2);
//////        bigimgb(rgy($:-1:1),rgx)=img(:,:,3);
//////        
//////    else
//////        //not a perfect match we need to interpolate
//////    end
//////    
//////end
//////
//////if atomsIsInstalled('SIVP') then
//////    imshow(bigimg);
//////end
//////RGBbigimg=hypermat([size(yy,2),size(xx,2),3],[bigimgr,bigimgg,bigimgb]);
//////
//////
//////if atomsIsInstalled('SIVP') then
//////    imshow(RGBbigimg);
//////end
//////
//////
//////
//////
//////h0=scf();
//////h0.color_map=bonecolormap(256);
//////Matplot(RGBbigimg(:,:,1));
//////
//////xc=140;yc=670;wc=1024*3+840,hc=768*3;
//////if atomsIsInstalled('SIVP') then
//////    imshow(imcrop(RGBbigimg,[xc,yc,wc,hc]));
//////end
//////
//////
////////olsfile1='/usr/local/diodesLASER/Techno/Plaques LAAS/2107/2107-28/ap-enduction-lext/2107-28-ap-enduction004_01a.ols';
////////olsfile2='/usr/local/diodesLASER/Techno/Plaques LAAS/2107/2107-25/ap-icp/LEXT/2107-25_006.ols';
////////olsfile3='/usr/local/diodesLASER/Techno/Plaques LAAS/2107/2107-25/ap-icp/LEXT/2107-25_008.ols';
////////olsfile4='/usr/local/diodesLASER/Techno/Plaques LAAS/2107/2107-25/ap-icp/LEXT/2107-25_040.ols';
//////
