diff --git a/@MappedTensor/arrayfun.m b/@MappedTensor/arrayfun.m index c263b06..181d65f 100644 --- a/@MappedTensor/arrayfun.m +++ b/@MappedTensor/arrayfun.m @@ -290,9 +290,13 @@ end % handle early return with single chunk - if bEarlyReturn && ~isempty(tData) && any(tData) - mtNewVar = tData; - return; + if bEarlyReturn + if ~isempty(tData) && any(tData(:)) + mtNewVar = tData; + return; + else + continue; + end end % - Write results diff --git a/@MappedTensor/arrayfun2.m b/@MappedTensor/arrayfun2.m index 5d610e1..3ef71d3 100644 --- a/@MappedTensor/arrayfun2.m +++ b/@MappedTensor/arrayfun2.m @@ -256,9 +256,13 @@ end % handle early return with single chunk - if bEarlyReturn && ~isempty(tData1) && any(tData1) - mtNewVar = tData1; - return; + if bEarlyReturn + if ~isempty(tData1) && any(tData1) + mtNewVar = tData1; + return; + else + continue; + end end % - Write results diff --git a/@MappedTensor/load.m b/@MappedTensor/load.m index a338d37..c5fdbc5 100644 --- a/@MappedTensor/load.m +++ b/@MappedTensor/load.m @@ -4,11 +4,14 @@ % | Extension | Description | % |-------------------|---------------------------| % | EDF | ESRF Data Format | -% | POS | Atom Probe Tomography | -% | NPY | Python NumPy array | +% | IMG SMV | ADSC X-ray detector image | +% | MAR MCCD | MAR CCD image | % | MRC MAP CCP4 RES | MRC MRC/CCP4/MAP electronic density map | -% | MAR | MAR CCD image | -% | IMG MCCD | ADSC X-ray detector image | +% | NRRD | Nearly Raw Raster Data | +% | NPY | Python NumPy array | +% | POS | Atom Probe Tomography | +% | VOL (PAR) | PyHST2 volume (tomography)| +% % % Not implemented % VOL + PAR PyHST2 volume reconstruction @@ -77,20 +80,26 @@ args = {}; header = []; switch upper(e) case '.EDF' - [Descr,args] = private_load_edf(filename); -case '.POS' - [Descr,args] = private_load_pos(filename); -case '.NPY' - [Descr,args,header] = private_load_npy(filename); -case {'.MRC','.MAP','.CCP4','.RES'} - [Descr,args,header] = private_load_mrc(filename); + [Descr,args,header] = private_load_edf(filename); +case '.IMG' + % ADSC (IMG/SMV) + [Descr,args,header] = private_load_smv(filename); + % could be a Mayo Clinic Analyze when supported... case {'.MAR','.MCCD'} [Descr,args,header] = private_load_mar(filename); +case {'.MRC','.MAP','.CCP4','.RES'} + [Descr,args,header] = private_load_mrc(filename); +case '.NPY' + [Descr,args,header] = private_load_npy(filename); +case '.NRRD' + [Descr,args,header] = private_load_nrrd(filename); +case '.POS' + [Descr,args,header] = private_load_pos(filename); case '.SMV' % ADSC (IMG/SMV) [Descr,args,header] = private_load_smv(filename); -case '.IMG' - % ADSC (IMG/SMV) - [Descr,args,header] = private_load_smv(filename); +case '.VOL' + [Descr,args,header] = private_load_vol(filename); + otherwise error([ mfilename ': unsupported file type ' upper(e) ' for ' filename ]); end diff --git a/@MappedTensor/plot.m b/@MappedTensor/plot.m index 9f70a47..8251047 100644 --- a/@MappedTensor/plot.m +++ b/@MappedTensor/plot.m @@ -3,6 +3,10 @@ % H = PLOT(M) plots the array M as a vector, matrix/surface, or volume. % Higher dimension tensors plot a projection. % +% PLOT(...,'log') transforms the signal in log-scale before plotting. +% +% PLOT(..., 'PROP1',VALUE1,...) uses specified plot properties. +% % Example: h=plot(MappedTensor(rand(100))); close(gcf); 1 h = []; @@ -22,6 +26,16 @@ end signal = squeeze(double(signal)); +% handle special 'log' option +for index=1:numel(varargin) + if ischar(varargin{index}) && strcmp(varargin{index}, 'log') + varargin(index) = []; + signal = signal - min(signal(:)); + signal(signal <=0) = nan; + signal = log10(signal); + end +end + if ~isempty(inputname(1)) iname = inputname(1); else @@ -39,9 +53,9 @@ % call the single plot method try - h = iFunc_plot(name, signal); + h = single_plot(name, signal, varargin{:}); if ~isempty(h) - h = iFunc_plot_menu(h, a, name); + h = single_plot_menu(h, a, name); end catch ME warning(getReport(ME)) @@ -51,7 +65,7 @@ % ------------------------------------------------------------------------------ % simple plot of the array -function h=iFunc_plot(name, signal, varargin) +function h=single_plot(name, signal, varargin) % this internal function plots a single model, 1D, 2D or 3D. h=[]; if isempty(signal) || isscalar(signal) @@ -74,7 +88,7 @@ set(h, 'DisplayName', name); %------------------------------------------------------------------------------- -function h=iFunc_plot_menu(h, a, name) +function h=single_plot_menu(h, a, name) % contextual menu for the single object being displayed % internal functions must be avoided as it uses LOTS of memory diff --git a/@MappedTensor/private/private_load_nrrd.m b/@MappedTensor/private/private_load_nrrd.m new file mode 100644 index 0000000..c3703db --- /dev/null +++ b/@MappedTensor/private/private_load_nrrd.m @@ -0,0 +1,214 @@ +function [Descr, args, header] = private_load_nrrd(filename) +% PRIVATE_LOAD_NRRD Read NRRD Nearly Raw Raster Data. +% +% A NRRD file is a simple binary block, preceeded by a text header until an +% empty line. +% +% See http://teem.sourceforge.net/nrrd/format.html + +Descr=''; args = []; + +header = nrrdread(filename); + +Descr= 'Nearly Raw Raster Data'; +args = { + 'Offset', header.Offset, ... + 'Format', header.Format, ... + 'MachineFormat', header.MachineFormat, ... + 'Dimension' , header.Dimension' }; + +% ------------------------------------------------------------------------------ +% https://www.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader +% ------------------------------------------------------------------------------ + +% Copyright (c) 2016, The MathWorks, Inc. +% All rights reserved. + +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: + +% * Redistributions of source code must retain the above copyright notice, this +% list of conditions and the following disclaimer. + +% * Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% * In all cases, the software is, and all modifications and derivatives of the +% software shall be, licensed to you solely for use in conjunction with +% MathWorks products and service offerings. + +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE +% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +function [meta] = nrrdread(filename) +%NRRDREAD Import NRRD imagery and metadata. +% [X, META] = NRRDREAD(FILENAME) reads the image volume and associated +% metadata from the NRRD-format file specified by FILENAME. +% +% Current limitations/caveats: +% * "Block" datatype is not supported. +% * Only tested with "gzip" and "raw" file encodings. +% * Very limited testing on actual files. +% * I only spent a couple minutes reading the NRRD spec. +% +% See the format specification online: +% http://teem.sourceforge.net/nrrd/format.html +% Copyright 2012 The MathWorks, Inc. +% Open file. +fid = fopen(filename, 'rb'); +assert(fid > 0, 'Could not open file.'); +cleaner = onCleanup(@() fclose(fid)); +% Magic line. +theLine = fgetl(fid); +assert(numel(theLine) >= 4, 'Bad signature in file.') +assert(isequal(theLine(1:4), 'NRRD'), 'Bad signature in file.') +% The general format of a NRRD file (with attached header) is: +% +% NRRD000X +% : +% : +% # +% ... +% : +% := +% := +% := +% # +% +% ... +meta = struct([]); +% Parse the file a line at a time. +while (true) + theLine = fgetl(fid); + + if (isempty(theLine) || feof(fid)) + % End of the header. + break; + end + + if (isequal(theLine(1), '#')) + % Comment line. + continue; + end + + % "fieldname:= value" or "fieldname: value" or "fieldname:value" + parsedLine = regexp(theLine, ':=?\s*', 'split','once'); + + assert(numel(parsedLine) == 2, 'Parsing error') + + field = lower(parsedLine{1}); + value = parsedLine{2}; + + field(isspace(field)) = ''; + meta(1).(field) = value; + +end +datatype = getDatatype(meta.type); +% Get the size of the data. +assert(isfield(meta, 'sizes') && ... + isfield(meta, 'encoding') , ... + 'Missing required metadata fields.') +dims = sscanf(meta.sizes, '%d'); +if isfield(meta,'dimension') + ndims = sscanf(meta.dimension, '%d'); + assert(numel(dims) == ndims); +end +meta.Offset = ftell(fid); +%data = readData(fid, meta, datatype); +%data = adjustEndian(data, meta); +% Reshape and get into MATLAB's order. +meta.Dimension = dims; +%X = reshape(data, dims'); +%X = permute(X, [2 1 3]); + +% prepare for MappedTensor +meta.Format = datatype; +meta.MachineFormat = ''; +if isfield(meta, 'endian') + switch meta.endian + case {'little','L','le'} + meta.MachineFormat = 'ieee-le'; + case {'big','B','be'} + meta.MachineFormat = 'ieee-be'; + end +end + +% ------------------------------------------------------------------------------ + +function datatype = getDatatype(metaType) +% Determine the datatype +switch (metaType) + case {'signed char', 'int8', 'int8_t'} + datatype = 'int8'; + + case {'uchar', 'unsigned char', 'uint8', 'uint8_t'} + datatype = 'uint8'; + case {'short', 'short int', 'signed short', 'signed short int', ... + 'int16', 'int16_t'} + datatype = 'int16'; + + case {'ushort', 'unsigned short', 'unsigned short int', 'uint16', ... + 'uint16_t'} + datatype = 'uint16'; + + case {'int', 'signed int', 'int32', 'int32_t'} + datatype = 'int32'; + + case {'uint', 'unsigned int', 'uint32', 'uint32_t'} + datatype = 'uint32'; + + case {'longlong', 'long long', 'long long int', 'signed long long', ... + 'signed long long int', 'int64', 'int64_t'} + datatype = 'int64'; + + case {'ulonglong', 'unsigned long long', 'unsigned long long int', ... + 'uint64', 'uint64_t'} + datatype = 'uint64'; + + case {'float'} + datatype = 'single'; + + case {'double'} + datatype = 'double'; + + otherwise + assert(false, 'Unknown datatype') +end + + +function data = readData(fidIn, meta, datatype) +switch (meta.encoding) + case {'raw'} + + data = fread(fidIn, inf, [datatype '=>' datatype]); + + case {'gzip', 'gz'} + error('Compressed NRRD is not supported') + + case {'txt', 'text', 'ascii'} + + data = fscanf(fidIn, '%f'); + data = cast(data, datatype); + + otherwise + assert(false, 'Unsupported encoding') +end + + +function data = adjustEndian(data, meta) +if ~isfield(meta, 'endian'), return; end +[~,~,endian] = computer(); +needToSwap = (isequal(endian, 'B') && isequal(lower(meta.endian), 'little')) || ... + (isequal(endian, 'L') && isequal(lower(meta.endian), 'big')); + +if (needToSwap) + data = swapbytes(data); +end diff --git a/@MappedTensor/private/private_load_vol.m b/@MappedTensor/private/private_load_vol.m new file mode 100644 index 0000000..697e568 --- /dev/null +++ b/@MappedTensor/private/private_load_vol.m @@ -0,0 +1,70 @@ +function [Descr, args, header] = private_load_vol(fname) +% PRIVATE_LOAD_VOL Read a RAW volume (float32, tomography/PyHST2)) +% +% Read a RAW volume (float32) generated by PyHST2: +% +% +% +% A '.par' file will be searched at the same location, to extract the volume shape. +% This '.par' file should contain text with lines such as: +% END_VOXEL_1 = 693 # X-end of reconstruction volume +% END_VOXEL_2 = 693 # Y-end of reconstruction volume +% END_VOXEL_3 = 961 # Z-end of reconstruction volume +% or +% NUM_IMAGE_1 = 693 # Number of pixels horizontally +% NUM_IMAGE_2 = 961 # Number of pixels vertically +% +% References: +% Alessandro Mirone et al, NIM B 324 (2014) 41. DOI: 10.1016/j.nimb.2013.09.030 + + +Descr=''; args = {}; frame = []; + +% get any PAR file as header +[p,f,e] = fileparts(fname); +shape = [ 0 0 0 ]; + +d = dir(fullfile(p, '*.par')); +if ~isempty(d) + par = fullfile(p, d(1).name); + disp([ mfilename ': Using PAR file ' par ]); + header = fileread(par); + header= str2struct(header); + % now search for fields and extract + + f = { 'END_VOXEL_1','END_VOXEL_2','END_VOXEL_3','NUM_IMAGE_1','NUM_IMAGE_1'}; + for index=1:numel(f) + field = f{index}; + if isfield(header, field) + if strcmp(field(end-1:end),'_1') dim=1; + elseif strcmp(field(end-1:end),'_2') dim=2; + elseif strcmp(field(end-1:end),'_3') dim=3; + end + if isnumeric(header.(field)) + val = header.(field); + else + val = str2num(strtok(header.(field),'#;%/')); + end + if isfinite(val) && ~shape(dim), shape(dim) = val; end + end + end + if numel(find(~shape)) == 1 + shape(~shape) = numel(frame.data)/prod(shape(shape ~= 0)); + end + header.shape = shape; +end + +header.Offset = 0; +header.Format = 'single'; % float32 +if all(shape > 0) + header.Dimension = shape; +else + % assume long vector + d = dir(fname); + header.Dimension = [ d.bytes/4 1 1 ]; +end + +Descr = 'RAW binary volume'; +args = { ... + 'Format', header.Format, ... + 'Dimension', header.Dimension }; diff --git a/README.md b/README.md index 99418b6..2dd97e6 100644 --- a/README.md +++ b/README.md @@ -171,11 +171,13 @@ with the following data formats. | Extension | Description | |-------------------|---------------------------| | EDF | ESRF Data Format (2D) | -| POS | Atom Probe Tomography (4 columns) | -| NPY | Python NumPy array (nD) | -| MRC MAP CCP4 RES | MRC MRC/CCP4/MAP electronic density map (3D) | -| MAR | MAR CCD image (2D) | -| IMG MCCD | ADSC X-ray detector image SMV (2D) | +| IMG SMV | [ADSC X-ray detector image SMV (2D)](https://strucbio.biologie.uni-konstanz.de/ccp4wiki/index.php/SMV_file_format) | +| MAR MCCD | MAR CCD image (2D) | +| MRC MAP CCP4 RES | [MRC MRC/CCP4/MAP electronic density map (3D)](https://www.ccpem.ac.uk/mrc_format/mrc2014.php) | +| NPY | [Python NumPy array (nD)](https://numpy.org/devdocs/reference/generated/numpy.lib.format.html) | +| NRRD | [Nearly Raw Raster Data](http://teem.sourceforge.net/nrrd/index.html) | +| POS | [Atom Probe Tomography (4 columns)](https://github.com/oscarbranson/apt-tools/blob/master/file-format-info.pdf) | +| VOL (PAR) | [PyHST2 tomography volume (3D)](http://ftp.esrf.fr/scisoft/PYHST2/output_format.html) | ## Using the array