From 2a612d2ac89fedf1d7728b124bd3c72b77f69c5e Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 18:25:03 +0200 Subject: [PATCH 1/6] load: add VOL import --- @MappedTensor/load.m | 3 + @MappedTensor/private/private_load_vol.m | 70 ++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 @MappedTensor/private/private_load_vol.m diff --git a/@MappedTensor/load.m b/@MappedTensor/load.m index a338d37..d95d92f 100644 --- a/@MappedTensor/load.m +++ b/@MappedTensor/load.m @@ -9,6 +9,7 @@ % | MRC MAP CCP4 RES | MRC MRC/CCP4/MAP electronic density map | % | MAR | MAR CCD image | % | IMG MCCD | ADSC X-ray detector image | +% | VOL | RAW volume (tomography) | % % Not implemented % VOL + PAR PyHST2 volume reconstruction @@ -88,6 +89,8 @@ [Descr,args,header] = private_load_mar(filename); case '.SMV' % ADSC (IMG/SMV) [Descr,args,header] = private_load_smv(filename); +case '.VOL' + [Descr,args,header] = private_load_vol(filename); case '.IMG' % ADSC (IMG/SMV) [Descr,args,header] = private_load_smv(filename); 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 }; From c513df5b7ff9698d143f7eb15700ef52e2fbe0fc Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 18:25:30 +0200 Subject: [PATCH 2/6] plot:add 'log' plot option. --- @MappedTensor/plot.m | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/@MappedTensor/plot.m b/@MappedTensor/plot.m index 9f70a47..f947b27 100644 --- a/@MappedTensor/plot.m +++ b/@MappedTensor/plot.m @@ -3,6 +3,8 @@ % H = PLOT(M) plots the array M as a vector, matrix/surface, or volume. % Higher dimension tensors plot a projection. % +% PLOT(..., 'PROP1',VALUE1,...) uses specified plot properties. +% % Example: h=plot(MappedTensor(rand(100))); close(gcf); 1 h = []; @@ -22,6 +24,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 +51,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 +63,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 +86,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 From a856ab1b24ad4676616a63f2c752faf2fe56007b Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 18:26:30 +0200 Subject: [PATCH 3/6] arrayfun: fix bug with EarlyReturn to continue along all chunck until 'true' --- @MappedTensor/arrayfun.m | 10 +++++++--- @MappedTensor/arrayfun2.m | 10 +++++++--- @MappedTensor/plot.m | 2 ++ 3 files changed, 16 insertions(+), 6 deletions(-) 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/plot.m b/@MappedTensor/plot.m index f947b27..8251047 100644 --- a/@MappedTensor/plot.m +++ b/@MappedTensor/plot.m @@ -3,6 +3,8 @@ % 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 From 4a925d18cdbf030de058aed9de385b0a9d1e5711 Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 19:19:57 +0200 Subject: [PATCH 4/6] update doc with load/VOL format --- @MappedTensor/load.m | 2 +- README.md | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/@MappedTensor/load.m b/@MappedTensor/load.m index d95d92f..c153bf2 100644 --- a/@MappedTensor/load.m +++ b/@MappedTensor/load.m @@ -9,7 +9,7 @@ % | MRC MAP CCP4 RES | MRC MRC/CCP4/MAP electronic density map | % | MAR | MAR CCD image | % | IMG MCCD | ADSC X-ray detector image | -% | VOL | RAW volume (tomography) | +% | VOL (PAR) | RAW volume (tomography) | % % Not implemented % VOL + PAR PyHST2 volume reconstruction diff --git a/README.md b/README.md index 99418b6..f8d1904 100644 --- a/README.md +++ b/README.md @@ -176,6 +176,7 @@ with the following data formats. | 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) | +| VOL (PAR) | PyHST2 tomography volume (3D) | ## Using the array From 9c2c9898f84b21f13dc775d50fa67ed5077b0f7b Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 21:40:08 +0200 Subject: [PATCH 5/6] Add NRRD lazy loader --- @MappedTensor/load.m | 36 ++-- @MappedTensor/private/private_load_nrrd.m | 230 ++++++++++++++++++++++ README.md | 13 +- 3 files changed, 258 insertions(+), 21 deletions(-) create mode 100644 @MappedTensor/private/private_load_nrrd.m diff --git a/@MappedTensor/load.m b/@MappedTensor/load.m index c153bf2..c5fdbc5 100644 --- a/@MappedTensor/load.m +++ b/@MappedTensor/load.m @@ -4,12 +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 | -% | VOL (PAR) | RAW volume (tomography) | +% | 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 @@ -78,22 +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 '.VOL' [Descr,args,header] = private_load_vol(filename); -case '.IMG' - % ADSC (IMG/SMV) - [Descr,args,header] = private_load_smv(filename); + otherwise error([ mfilename ': unsupported file type ' upper(e) ' for ' filename ]); end diff --git a/@MappedTensor/private/private_load_nrrd.m b/@MappedTensor/private/private_load_nrrd.m new file mode 100644 index 0000000..e3412aa --- /dev/null +++ b/@MappedTensor/private/private_load_nrrd.m @@ -0,0 +1,230 @@ +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'} + tmpBase = tempname(); + tmpFile = [tmpBase '.gz']; + fidTmp = fopen(tmpFile, 'wb'); + assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression') + + tmp = fread(fidIn, inf, 'uint8=>uint8'); + fwrite(fidTmp, tmp, 'uint8'); + fclose(fidTmp); + + gunzip(tmpFile) + + fidTmp = fopen(tmpBase, 'rb'); + cleaner = onCleanup(@() fclose(fidTmp)); + + meta.encoding = 'raw'; + data = readData(fidTmp, meta, datatype); + + 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/README.md b/README.md index f8d1904..2dd97e6 100644 --- a/README.md +++ b/README.md @@ -171,12 +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) | -| VOL (PAR) | PyHST2 tomography volume (3D) | +| 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 From 41dddbefa77b53c043e980d71eefdcf237870008 Mon Sep 17 00:00:00 2001 From: Emmanuel FARHI Date: Mon, 18 May 2020 21:46:39 +0200 Subject: [PATCH 6/6] load: NRRD: gzip is unsupported. --- @MappedTensor/private/private_load_nrrd.m | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/@MappedTensor/private/private_load_nrrd.m b/@MappedTensor/private/private_load_nrrd.m index e3412aa..c3703db 100644 --- a/@MappedTensor/private/private_load_nrrd.m +++ b/@MappedTensor/private/private_load_nrrd.m @@ -10,7 +10,6 @@ header = nrrdread(filename); - Descr= 'Nearly Raw Raster Data'; args = { 'Offset', header.Offset, ... @@ -192,22 +191,7 @@ data = fread(fidIn, inf, [datatype '=>' datatype]); case {'gzip', 'gz'} - tmpBase = tempname(); - tmpFile = [tmpBase '.gz']; - fidTmp = fopen(tmpFile, 'wb'); - assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression') - - tmp = fread(fidIn, inf, 'uint8=>uint8'); - fwrite(fidTmp, tmp, 'uint8'); - fclose(fidTmp); - - gunzip(tmpFile) - - fidTmp = fopen(tmpBase, 'rb'); - cleaner = onCleanup(@() fclose(fidTmp)); - - meta.encoding = 'raw'; - data = readData(fidTmp, meta, datatype); + error('Compressed NRRD is not supported') case {'txt', 'text', 'ascii'}