[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: How to use debugging tools
|
From: |
Colin Ingram |
|
Subject: |
Re: How to use debugging tools |
|
Date: |
Mon, 19 Sep 2005 19:56:23 -0500 |
|
User-agent: |
Mozilla Thunderbird 1.0 (Windows/20041206) |
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
-------------------------------------------------------------
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