From help-request at octave dot org Mon Sep 19 19:57:05 2005 Subject: Re: How to use debugging tools From: Colin Ingram To: help at octave dot org Date: Mon, 19 Sep 2005 19:56:23 -0500 This is a multi-part message in MIME format. --------------020609090105050903070303 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Colin Ingram wrote: > > ************************************************************** > I've attached tiffread.m. I've experienced some more strange behavior > with the db tools but we will save that for later. > > Oh bother I've really attached it this time. > > ------------------------------------------------------------- > Octave is freely available under the terms of the GNU GPL. > > Octave's home on the web: http://www.octave.org > How to fund new projects: http://www.octave.org/funding.html > Subscription information: http://www.octave.org/archive.html > ------------------------------------------------------------- > --------------020609090105050903070303 Content-Type: text/plain; name="tiffread.m" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="tiffread.m" function [stack, img_read] = tiffread(filename, img_first, img_last) % tiffread, version 2.1 % % [stack, nbImages] = tiffread; % [stack, nbImages] = tiffread(filename); % [stack, nbImages] = tiffread(filename, imageIndex); % [stack, nbImages] = tiffread(filename, firstImageIndex, lastImageIndex); % % Reads 8,16,32 bits uncompressed grayscale and (some) color tiff files, % as well as stacks or multiple tiff images, for example those produced % by metamorph or NIH-image. However, the entire TIFF standard is not % supported (but you may extend it). % % The function can be called with a file name in the current directory, % or without argument, in which case it pop up a file openning dialog % to allow manual selection of the file. % If the stacks contains multiples images, loading can be restricted by % specifying the first and last images to read, or just one image to read. % % at return, nbimages contains the number of images read, and S is a vector % containing the different images with some additional informations. The % image pixels values are stored in the field .data, for gray level images, % or in the fields .red, .green and .blue % the pixels values are in the native (integer) format, % and must be converted to be used in most matlab functions. % % Example: % im = tiffread('spindle.stk'); % imshow( double(im(5).data), [] ); % % Francois Nedelec, EMBL, Copyright 1999-2004. % Last modified July 7th, 2004 at Woods Hole during the physiology course. % Thanks to Kendra Burbank for suggesting the waitbar % % Please, help us improve this software: send us feedback/bugs/suggestions % This software is provided at no cost by a public research institution. % However, postcards are always welcome! % % Francois Nedelec % Cell Biology and Biophysics, EMBL; Meyerhofstrasse 1; 69117 Heidelberg; Germany % http://www.embl.org % http://www.cytosim.org ##Optimization: join adjacent TIF strips: this results in faster reads consolidateStrips = true; ##if there is no argument, we ask the user to choose a file: if (nargin == 0) ## **to-fix** [CJI: no octave ui compat.] ## [filename, pathname] = uigetfile('*.tif;*.stk;*.lsm', 'select image file'); ## filename = [ pathname, filename ]; fflush(stdout); filename=input("Input Filename?", "s") endif if (nargin<=1) img_first = 1; img_last = 10000; endif if (nargin==2) img_last = img_first; endif % not all valid tiff tags have been included, as they are really a lot... % if needed, tags can easily be added to this code % See the official list of tags: % http://partners.adobe.com/asn/developer/pdfs/tn/TIFF6.pdf % % the structure IMG is returned to the user, while TIF is not. % so tags usefull to the user should be stored as fields in IMG, while % those used only internally can be stored in TIF. global TIF; TIF = []; ## counters for the number of images read and skipped img_skip = 0; img_read = 0; ## set defaults values : TIF.SampleFormat = 1; TIF.SamplesPerPixel = 1; TIF.BOS ="l"; # byte order string if isempty(findstr(filename,".")) filename = [filename,".stk"]; endif [TIF.file, message] = fopen(filename,"rb","ieee-le"); if TIF.file == -1 filename = strrep(filename, ".stk", ".tif"); [TIF.file, message] = fopen(filename,"rb","ieee-le"); if TIF.file == -1 error(["Could not open <",filename,"> for reading!"]); endif endif ## 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 = "ieee-le"; # normal PC format elseif ( strcmp(byte_order',"MM") ) TIF.BOS = "ieee-be"; else error("This is not a TIFF file. Headers Specifies an Invalid File Format."); endif ##----- read in a number which identifies file as TIFF format tiff_id = fread(TIF.file,1,"uint16", 0, TIF.BOS); if (tiff_id ~= 42) error("This is not a TIFF file. Header Specifies an invalid TIFF file_ID"); endif ##----- read the byte offset for the first image file directory (IFD) ifd_pos = fread(TIF.file,1,"uint32", 0 ,TIF.BOS); while (ifd_pos ~= 0) clear IMG; IMG.filename = filename; ## move in the file to the first IFD fseek(TIF.file, ifd_pos, -1); ## read in the number of IFD entries num_entries = fread(TIF.file,1,"uint16", 0, TIF.BOS); ## 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", 0, TIF.BOS); entry = readIFDentry; switch TIF.entry_tag case 254 TIF.NewSubfiletype = entry.val; case 256 # image width - number of column IMG.width = entry.val; case 257 # image height - number of row IMG.height = entry.val; TIF.ImageLength = entry.val; case 258 # BitsPerSample per sample TIF.BitsPerSample = entry.val; TIF.BytesPerSample = TIF.BitsPerSample / 8; IMG.bits = TIF.BitsPerSample(1); case 259 # compression if (entry.val ~= 1) error("Compression format not supported."); endif case 262 % photometric interpretation TIF.PhotometricInterpretation = entry.val; if ( TIF.PhotometricInterpretation == 3 ) fprintf(1, "warning: ignoring the look-up table defined in the TIFF file"); endif case 269 IMG.document_name = entry.val; case 270 % comment: % TIF.info = entry.val; TIF.info = parseMM_info(entry.val); case 271 IMG.make = entry.val; case 273 % strip offset TIF.StripOffsets = entry.val; TIF.StripNumber = entry.cnt; case 277 % sample_per pixel TIF.SamplesPerPixel = entry.val; 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 282 % X resolution IMG.x_resolution = entry.val; case 283 % Y resolution IMG.y_resolution = 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 ) error(sprintf("PlanarConfiguration = %i not supported", TIF.PlanarConfiguration)); endif case 296 % resolution unit IMG.resolution_unit= entry.val; case 305 % software IMG.software = entry.val; case 306 % datetime IMG.datetime = entry.val; case 315 IMG.artist = entry.val; case 317 %predictor for compression if (entry.val ~= 1) error("unsuported predictor value"); end case 320 % color map IMG.cmap = entry.val; IMG.colors = entry.cnt/3; case 339 TIF.SampleFormat = entry.val; if ( TIF.SampleFormat > 2 ) error(sprintf("unsuported sample format = %i", TIF.SampleFormat)); endif case 33628 %metamorph specific data IMG.MM_private1 = entry.val; case 33629 %this tag identify the image as a Metamorph stack! TIF.MM_stack = entry.val; TIF.MM_stackCnt = entry.cnt; ## **to fix** GUI ## if ( img_last > img_first ) ## waitbar_handle = waitbar(0,"Please wait...","Name",["Reading " filename]); ## endif case 33630 %metamorph stack data: wavelength TIF.MM_wavelength = entry.val; case 33631 %metamorph stack data: gain/background? TIF.MM_private2 = entry.val; case 34412 % Zeiss LSM data (I have no idea what that represents...) IMG.LSM = entry.val; otherwise fprintf(1,"ignored TIFF entry with tag %i (cnt %i)\n", TIF.entry_tag, entry.cnt); endswitch # IMG = assignIFDentry(entry); # keyboard; ## move to next IFD entry in the file fseek(TIF.file, file_pos+12,-1); endfor ## total number of bytes per image: PlaneBytesCnt = IMG.width * IMG.height * TIF.BytesPerSample; 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; endif endwhile ## Consolidate the Strips if ( BytesCnt <= PlaneBytesCnt(1) ) && ( ConsolidateCnt > 1 ) TIF.StripByteCounts = [BytesCnt; TIF.StripByteCounts(ConsolidateCnt+1:TIF.StripNumber) ]; TIF.StripOffsets = TIF.StripOffsets( [1 , ConsolidateCnt+1:TIF.StripNumber]); TIF.StripNumber = 1 + TIF.StripNumber - ConsolidateCnt; endif endif endif %read the next IFD address: ifd_pos = fread(TIF.file, 1, "uint32", 0,TIF.BOS); %if (ifd_pos) disp(["next ifd at", num2str(ifd_pos)]); end if isfield( TIF, "MM_stack" ) if ( img_last > TIF.MM_stackCnt ) img_last = TIF.MM_stackCnt; end %this loop is to read metamorph stacks: 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 ) IMG.data = 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); 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); # if exist('waitbar_handle', 'var') # waitbar( img_read/TIF.MM_stackCnt, waitbar_handle); # end [ IMG.info, IMG.MM_stack, IMG.MM_wavelength, IMG.MM_private2 ] = extractMetamorphData(ii); img_read = img_read + 1; 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 ) IMG.data = read_plane(0, IMG.width, IMG.height, 1); 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); end img_read = img_read + 1; try stack( img_read ) = IMG; catch stack IMG error("The file contains dissimilar images: you can only read them one by one"); end else img_skip = img_skip + 1; end if ( img_skip + img_read >= img_last ) break; end end end %clean-up fclose(TIF.file); #if exist("waitbar_handle", "var") # delete( waitbar_handle ); # clear waitbar_handle; #end %return empty array if nothing was read if ~ exist( "stack", "var") stack = []; 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; endif %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 ); 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 % transpose the image plane = 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, 0,TIF.BOS ); if ( length(bytes) ~= StripLength ) error("End of file reached unexpectedly."); end strip = reshape(bytes, width, StripLength / width); return; %===================sub-functions that reads an IFD entry:=================== function [nbbytes, typechar] = octave_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"; case 6 nbbytes=1; typechar="int8"; case 8 nbbytes=2; typechar="int16"; case 9 nbbytes=4; typechar="int32"; case 10 nbbytes=8; typechar="int32"; case 11 nbbytes=4; typechar="single"; case 12 nbbytes=8; typechar="double"; otherwise error("tiff type not supported") endswitch return; endfunction %===================sub-functions that reads an IFD entry:=================== function entry = readIFDentry() global TIF; entry.typecode = fread(TIF.file, 1, "uint16", 0,TIF.BOS); entry.cnt = fread(TIF.file, 1, "uint32", 0,TIF.BOS); [ entry.nbbytes, entry.typechar ] = octave_type(entry.typecode); if entry.nbbytes * entry.cnt > 4 ## next field contains an offset: offset = fread(TIF.file, 1, "uint32", 0,TIF.BOS); fseek(TIF.file, offset, -1); endif if TIF.entry_tag == 33629 # special metamorph "rationals" entry.val = fread(TIF.file, 6*entry.cnt, entry.typechar, ... 0,TIF.BOS); elseif (entry.typechar == 5 || entry.typechar == 10) entry.val = fread(TIF.file, 2*entry.cnt, entry.typechar, 0, ... TIF.BOS); else entry.val = fread(TIF.file, entry.cnt, entry.typechar, ... 0,TIF.BOS); endif if ( entry.typecode == 2 ) entry.val = char(entry.val'); endif return; endfunction ##==================sub-functions that assigns an IFD entry:=================== function IMG = assignIFDentry(entry) global TIF; switch TIF.entry_tag case 254 TIF.NewSubfiletype = entry.val; case 256 # image width - number of column IMG.width = entry.val; case 257 # image height - number of row IMG.height = entry.val; TIF.ImageLength = entry.val; case 258 # BitsPerSample per sample TIF.BitsPerSample = entry.val; TIF.BytesPerSample = TIF.BitsPerSample / 8; IMG.bits = TIF.BitsPerSample(1); case 259 # compression if (entry.val ~= 1) error("Compression format not supported."); endif case 262 % photometric interpretation TIF.PhotometricInterpretation = entry.val; if ( TIF.PhotometricInterpretation == 3 ) fprintf(1, "warning: ignoring the look-up table defined in the TIFF file"); endif case 269 IMG.document_name = entry.val; case 270 % comment: % TIF.info = entry.val; TIF.info = parseMM_info(entry.val); case 271 IMG.make = entry.val; case 273 % strip offset TIF.StripOffsets = entry.val; TIF.StripNumber = entry.cnt; case 277 % sample_per pixel TIF.SamplesPerPixel = entry.val; 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 282 % X resolution IMG.x_resolution = entry.val; case 283 % Y resolution IMG.y_resolution = 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 ) error(sprintf("PlanarConfiguration = %i not supported", TIF.PlanarConfiguration)); endif case 296 % resolution unit IMG.resolution_unit= entry.val; case 305 % software IMG.software = entry.val; case 306 % datetime IMG.datetime = entry.val; case 315 IMG.artist = entry.val; case 317 %predictor for compression if (entry.val ~= 1) error("unsuported predictor value"); end case 320 % color map IMG.cmap = entry.val; IMG.colors = entry.cnt/3; case 339 TIF.SampleFormat = entry.val; if ( TIF.SampleFormat > 2 ) error(sprintf("unsuported sample format = %i", TIF.SampleFormat)); endif case 33628 %metamorph specific data IMG.MM_private1 = entry.val; case 33629 %this tag identify the image as a Metamorph stack! TIF.MM_stack = entry.val; TIF.MM_stackCnt = entry.cnt; ## **to fix** GUI ## if ( img_last > img_first ) ## waitbar_handle = waitbar(0,"Please wait...","Name",["Reading " filename]); ## endif case 33630 %metamorph stack data: wavelength TIF.MM_wavelength = entry.val; case 33631 %metamorph stack data: gain/background? TIF.MM_private2 = entry.val; case 34412 % Zeiss LSM data (I have no idea what that represents...) IMG.LSM = entry.val; otherwise fprintf(1,"ignored TIFF entry with tag %i (cnt %i)\n", TIF.entry_tag, entry.cnt); endswitch %==============distribute the metamorph infos to each frame: function [info, stack, wavelength, private2 ] = extractMetamorphData(imgCnt) global TIF; info = []; stack = []; wavelength = []; private2 = []; if TIF.MM_stackCnt == 1 return; end left = imgCnt - 1; if isfield( TIF, "info" ) % S = length(TIF.info) / TIF.MM_stackCnt; info = TIF.info(imgCnt); end if isfield( TIF, "MM_stack" ) S = length(TIF.MM_stack) / TIF.MM_stackCnt; stack = TIF.MM_stack( [S*left+1:S*left+S] ); end if isfield( TIF, "MM_wavelength" ) S = length(TIF.MM_wavelength) / TIF.MM_stackCnt; wavelength = TIF.MM_wavelength( [S*left+1:S*left+S] ); end if isfield( TIF, "MM_private2" ) S = length(TIF.MM_private2) / TIF.MM_stackCnt; private2 = TIF.MM_private2( [S*left+1:S*left+S] ); end return; %%%% Parse the Metamorph camera info tag into respective fields % EVBR 2/7/2005 function mm_infor = parseMM_info(info) i=0; kk=info; kk(find(double(kk)==0)) = ""; while(length(kk)) i=i+1; [tmp, kk] = strtok(kk, 13); % [tmp,b] = strtok(tmp, 0); token{i} = sscanf(tmp, "%s"); end tok = {token{1:end-1}}; len = length(tok)/11; k=0; for i=1:length(tok) if ~isempty(token{i}) [tkk, remk] = strtok(token{i}, 58); if strcmp(tkk, "Exposure") k=k+1; end switch tkk case "Exposure" % Exposure: 200 ms -> 200 mm_infor(k).Exposure = str2double(remk(2:strfind(remk, "ms")-1)); case "Binning" % Binning: 1 x 1 -> [1 1] tmp = sscanf(remk, "%c %d %c %d")'; mm_infor(k).Binning = tmp([2,4]); case "Region" % Region: 1392 x 1040, offset at (0, 0) -> [1392 1040], [0 0] [a,b] = strtok(remk, ","); tmp = sscanf(a, "%c %d %c %d")'; mm_infor(k).Region.Size = tmp([2,4]); [a,b] = strtok(b, ","); mm_infor(k).Region.Offset = [str2double(a(end)), str2double(b(2))]; case "Subtract" % Subtract: Off -> 0 mm_infor(k).Subtract = ~strcmp(remk(2:end),"Off"); case "Shading" % Shading: Off -> 0 mm_infor(k).Shading = ~strcmp(remk(2:end),"Off"); case "Digitizer" % Digitizer: 5MHz -> 5 mm_infor(k).Digitizer = sscanf(remk(2:end), "%d %*s")'; case "Gain" % Gain: Gain 1 (1x) -> 1 mm_infor(k).Gain = sscanf(remk(strfind(remk,"(")+1:end),... "%d %*s %d")'; case "CameraShutter" % Camera Shutter: Always Open -> 'AlwaysOpen' mm_infor(k).CameraShutter = remk(2:end); case "ClearCount" % Clear Count: 2 -> 2 tmp = sscanf(remk, "%c %d")'; mm_infor(k).ClearCount = tmp(2); case "TriggerMode" % Trigger Mode: Normal -> "Normal" mm_infor(k).TriggerMode = remk(2:end); case "Temperature" % Temperature: 20 -> 20 tmp = sscanf(remk, "%c %d")'; mm_infor(k).Temperature = tmp(2); otherwise warning("Uknown MM_Tag") end end end --------------020609090105050903070303-- ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------