subsasgn

 Subscript assignment
 See subsref for meaning of fields.
 _______________________________________________________________________
 Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

0001 function obj = subsasgn(obj,subs,varargin)
0002 % Subscript assignment
0003 % See subsref for meaning of fields.
0004 % _______________________________________________________________________
0005 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0006 
0007 %
0008 % $Id: subsasgn.m 253 2005-10-13 15:31:34Z guillaume $
0009 
0010 
0011 switch subs(1).type,
0012 case {'.'},
0013     if numel(obj)~=nargin-2,
0014         error('The number of outputs should match the number of inputs.');
0015     end;
0016     objs = struct(obj);
0017     for i=1:length(varargin),
0018         val     = varargin{i};
0019         obji    = class(objs(i),'read_nifti');
0020         obji    = fun(obji,subs,val);
0021         objs(i) = struct(obji);
0022     end;
0023     obj = class(objs,'read_nifti');
0024 
0025 case {'()'},
0026     objs = struct(obj);
0027     if length(subs)>1,
0028         t    = subsref(objs,subs(1));
0029         % A lot of this stuff is a little flakey, and may cause Matlab to bomb.
0030         %
0031         %if numel(t) ~= nargin-2,
0032         %    error('The number of outputs should match the number of inputs.');
0033         %end;
0034         for i=1:numel(t),
0035             val  = varargin{1};
0036             obji = class(t(i),'read_nifti');
0037             obji = subsasgn(obji,subs(2:end),val);
0038             t(i) = struct(obji);
0039         end;
0040         objs = subsasgn(objs,subs(1),t);
0041     else
0042         if numel(varargin)>1,
0043             error('Illegal right hand side in assignment. Too many elements.');
0044         end;
0045         val = varargin{1};
0046         if isa(val,'read_nifti'),
0047             objs = subsasgn(objs,subs,struct(val));
0048         elseif isempty(val),
0049             objs = subsasgn(objs,subs,[]);
0050         else
0051             error('Assignment between unlike types is not allowed.');
0052         end;
0053     end;
0054     obj = class(objs,'read_nifti');
0055 
0056 otherwise
0057     error('Cell contents reference from a non-cell array object.');
0058 end;
0059 return;
0060 %=======================================================================
0061 
0062 %=======================================================================
0063 function obj = fun(obj,subs,val)
0064 % Subscript referencing
0065 
0066 switch subs(1).type,
0067 case {'.'},
0068     objs = struct(obj);
0069     for ii=1:numel(objs)
0070         obj = objs(ii);
0071 
0072         if any(strcmpi(subs(1).subs,{'dat'})),
0073             if length(subs)>1,
0074                 val = subsasgn(obj.dat,subs(2:end),val);
0075             end;
0076             obj      = assigndat(obj,val);
0077             objs(ii) = obj;
0078             continue;
0079         end;
0080 
0081         if isempty(obj.hdr), obj.hdr = empty_hdr; end;
0082         if ~isfield(obj.hdr,'magic'), error('Not a NIFTI-1 header'); end;
0083 
0084         if length(subs)>1, % && ~strcmpi(subs(1).subs,{'raw','dat'}),
0085             val0 = subsref(class(obj,'read_nifti'),subs(1));
0086             val1 = subsasgn(val0,subs(2:end),val);
0087         else
0088             val1 = val;
0089         end;
0090 
0091         switch(subs(1).subs)
0092         case {'extras'}
0093             if length(subs)>1,
0094                 obj.extras = subsasgn(obj.extras,subs(2:end),val);
0095             else
0096                 obj.extras = val;
0097             end;
0098 
0099         case {'mat0'}
0100             if ~isnumeric(val1) || ndims(val1)~=2 || any(size(val1)~=[4 4]) || sum((val1(4,:)-[0 0 0 1]).^2)>1e-8,
0101                 error('"mat0" should be a 4x4 matrix, with a last row of 0,0,0,1.');
0102             end;
0103             if obj.hdr.qform_code==0, obj.hdr.qform_code=2; end;
0104             s = double(bitand(obj.hdr.xyzt_units,7));
0105             if s
0106                 d = findindict(s,'units');
0107                 val1 = diag([[1 1 1]/d.rescale 1])*val1;
0108             end;
0109             obj.hdr = encode_qform0(double(val1), obj.hdr);
0110 
0111         case {'mat0_intent'}
0112             if isempty(val1),
0113                 obj.hdr.qform_code = 0;
0114             else
0115                 if ~ischar(val1) && ~(isnumeric(val1) && numel(val1)==1),
0116                     error('"mat0_intent" should be a string or a scalar.');
0117                 end;
0118                 d = findindict(val1,'xform');
0119                 if ~isempty(d)
0120                     obj.hdr.qform_code = d.code;
0121                 end;
0122             end;
0123 
0124         case {'mat'}
0125             if ~isnumeric(val1) || ndims(val1)~=2 || any(size(val1)~=[4 4]) || sum((val1(4,:)-[0 0 0 1]).^2)>1e-8
0126                 error('"mat" should be a 4x4 matrix, with a last row of 0,0,0,1.');
0127             end;
0128             if obj.hdr.sform_code==0, obj.hdr.sform_code=2; end;
0129             s = double(bitand(obj.hdr.xyzt_units,7));
0130             if s
0131                 d = findindict(s,'units');
0132                 val1 = diag([[1 1 1]/d.rescale 1])*val1;
0133             end;
0134             val1           = val1 * [eye(4,3) [1 1 1 1]'];
0135             obj.hdr.srow_x = val1(1,:);
0136             obj.hdr.srow_y = val1(2,:);
0137             obj.hdr.srow_z = val1(3,:);
0138 
0139         case {'mat_intent'}
0140             if isempty(val1),
0141                 obj.hdr.sform_code = 0;
0142             else
0143                 if ~ischar(val1) && ~(isnumeric(val1) && numel(val1)==1),
0144                     error('"mat_intent" should be a string or a scalar.');
0145                 end;
0146                 d = findindict(val1,'xform');
0147                 if ~isempty(d),
0148                     obj.hdr.sform_code = d.code;
0149                 end;
0150             end;
0151 
0152         case {'intent'}
0153             if ~valid_fields(val1,{'code','param','name'})
0154                 obj.hdr.intent_code = 0;
0155                 obj.hdr.intent_p1   = 0;
0156                 obj.hdr.intent_p2   = 0;
0157                 obj.hdr.intent_p3   = 0;
0158                 obj.hdr.intent_name = '';
0159             else
0160                 if ~isfield(val1,'code'),
0161                     val1.code = obj.hdr.intent_code;
0162                 end;
0163                 d = findindict(val1.code,'intent');
0164                 if ~isempty(d),
0165                     obj.hdr.intent_code = d.code;
0166                     if isfield(val1,'param'),
0167                         prm = [double(val1.param(:))  ; 0 ; 0; 0];
0168                         prm = [prm(1:length(d.param)) ; 0 ; 0; 0];
0169                         obj.hdr.intent_p1 = prm(1);
0170                         obj.hdr.intent_p2 = prm(2);
0171                         obj.hdr.intent_p3 = prm(3);
0172                     end;
0173                     if isfield(val1,'name'),
0174                         obj.hdr.intent_name = val1.name;
0175                     end;
0176                 end;
0177             end;
0178 
0179          case {'diminfo'}
0180             if ~valid_fields(val1,{'frequency','phase','slice','slice_time'})
0181                 tmp = obj.hdr.dim_info;
0182                 for bit=1:6,
0183                     tmp = bitset(tmp,bit,0);
0184                 end;
0185                 obj.hdr.dim_info       = tmp;
0186                 obj.hdr.slice_start    = 0;
0187                 obj.hdr.slice_end      = 0;
0188                 obj.hdr.slice_duration = 0;
0189                 obj.hdr.slice_code     = 0;
0190             else
0191                 if isfield(val1,'frequency'),
0192                     tmp = val1.frequency;
0193                     if ~isnumeric(tmp) || numel(tmp)~=1 || tmp<0 || tmp>3,
0194                         error('Invalid frequency direction');
0195                     end;
0196                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,1,bitget(tmp,1));
0197                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,2,bitget(tmp,2));
0198                 end;
0199 
0200                 if isfield(val1,'phase'),
0201                     tmp = val1.phase;
0202                     if ~isnumeric(tmp) || numel(tmp)~=1 || tmp<0 || tmp>3,
0203                         error('Invalid phase direction');
0204                     end;
0205                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,3,bitget(tmp,1));
0206                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,4,bitget(tmp,2));
0207                 end;
0208 
0209                 if isfield(val1,'slice'),
0210                     tmp = val1.slice;
0211                     if ~isnumeric(tmp) || numel(tmp)~=1 || tmp<0 || tmp>3,
0212                         error('Invalid slice direction');
0213                     end;
0214                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,5,bitget(tmp,1));
0215                     obj.hdr.dim_info = bitset(obj.hdr.dim_info,6,bitget(tmp,2));
0216                 end;
0217 
0218                 if isfield(val1,'slice_time')
0219                     tim = val1.slice_time;
0220                     if ~valid_fields(tim,{'start','end','duration','code'}),
0221                         obj.hdr.slice_code     = 0;
0222                         obj.hdr.slice_start    = 0;
0223                         obj.hdr.end_slice      = 0;
0224                         obj.hdr.slice_duration = 0;
0225                     else
0226                         % sld = double(bitget(obj.hdr.dim_info,5)) + 2*double(bitget(obj.hdr.dim_info,6));
0227 
0228                         if isfield(tim,'start'),
0229                             ss = double(tim.start);
0230                             if isnumeric(ss) && numel(ss)==1 && ~rem(ss,1), % && ss>=1 && ss<=obj.hdr.dim(sld+1)
0231                                 obj.hdr.slice_start = ss-1;
0232                             else
0233                                 error('Inappropriate "slice_time.start".');
0234                             end;
0235                         end;
0236 
0237                         if isfield(tim,'end'),
0238                             ss = double(tim.end);
0239                             if isnumeric(ss) && numel(ss)==1 && ~rem(ss,1), % && ss>=1 && ss<=obj.hdr.dim(sld+1)
0240                                 obj.hdr.slice_end = ss-1;
0241                             else
0242                                 error('Inappropriate "slice_time.end".');
0243                             end;
0244                         end;
0245 
0246                         if isfield(tim,'duration')
0247                             sd = double(tim.duration);
0248                             if isnumeric(sd) && numel(sd)==1,
0249                                 s  = double(bitand(obj.hdr.xyzt_units,24));
0250                                 d  = findindict(s,'units');
0251                                 if ~isempty(d) && d.rescale, sd = sd/d.rescale; end;
0252                                 obj.hdr.slice_duration = sd;
0253                             else
0254                                 error('Inappropriate "slice_time.duration".');
0255                             end;
0256                         end;
0257 
0258                         if isfield(tim,'code'),
0259                             d = findindict(tim.code,'sliceorder');
0260                             if ~isempty(d),
0261                                 obj.hdr.slice_code = d.code;
0262                             end;
0263                         end;
0264                     end;
0265                 end;
0266             end;
0267 
0268         case {'timing'}
0269             if ~valid_fields(val1,{'toffset','tspace'}),
0270                 obj.hdr.pixdim(5) = 0;
0271                 obj.hdr.toffset   = 0;
0272             else
0273                 s  = double(bitand(obj.hdr.xyzt_units,24));
0274                 d  = findindict(s,'units');
0275                 if isfield(val1,'toffset'),
0276                     if isnumeric(val1.toffset) && numel(val1.toffset)==1,
0277                         if d.rescale,
0278                             val1.toffset = val1.toffset/d.rescale;
0279                         end;
0280                         obj.hdr.toffset   = val1.toffset;
0281                     else
0282                         error('"timing.toffset" needs to be numeric with 1 element');
0283                     end;
0284                 end;
0285                 if isfield(val1,'tspace'),
0286                     if isnumeric(val1.tspace) && numel(val1.tspace)==1,
0287                         if d.rescale,
0288                             val1.tspace = val1.tspace/d.rescale;
0289                         end;
0290                         obj.hdr.pixdim(5) = val1.tspace;
0291                     else
0292                         error('"timing.tspace" needs to be numeric with 1 element');
0293                     end;
0294                 end;
0295             end;
0296 
0297         case {'descrip'}
0298             if isempty(val1), val1 = char(val1); end;
0299             if ischar(val1),
0300                 obj.hdr.descrip = val1;
0301             else
0302                 error('"descrip" must be a string.');
0303             end;
0304 
0305         case {'cal'}
0306             if isempty(val1),
0307                 obj.hdr.cal_min = 0;
0308                 obj.hdr.cal_max = 0;
0309             else
0310                 if isnumeric(val1) && numel(val1)==2,
0311                     obj.hdr.cal_min = val1(1);
0312                     obj.hdr.cal_max = val1(2);
0313                 else
0314                     error('"cal" should contain two elements.');
0315                 end;
0316             end;
0317 
0318         case {'aux_file'}
0319             if isempty(val1), val1 = char(val1); end;
0320             if ischar(val1),
0321                 obj.hdr.aux_file = val1;
0322             else
0323                 error('"aux_file" must be a string.');
0324             end;
0325 
0326         case {'hdr'}
0327             error('hdr is a read-only field.');
0328             obj.hdr = val1;
0329 
0330         otherwise
0331             error(['Reference to non-existent field ''' subs(1).subs '''.']);
0332         end;
0333 
0334         objs(ii) = obj;
0335     end
0336     obj = class(objs,'read_nifti');
0337 
0338 otherwise
0339     error('This should not happen.');
0340 end;
0341 return;
0342 %=======================================================================
0343 
0344 %=======================================================================
0345 function obj = assigndat(obj,val)
0346 if isa(val,'create_file_array'),
0347     sz = size(val);
0348     if numel(sz)>8,
0349         error('Too many dimensions in data.');
0350     end;
0351     sz = [sz 1 1 1 1 1 1 1 1];
0352     sz = sz(1:8);
0353     sval = struct(val);
0354     d    = findindict(sval.dtype,'dtype');
0355     if isempty(d)
0356         error(['Unknown datatype (' num2str(double(sval.datatype)) ').']);
0357     end;
0358 
0359     [pth,nam,suf]    = fileparts(sval.fname);
0360     if any(strcmp(suf,{'.img','.IMG'}))
0361         val.offset = max(sval.offset,0);
0362         obj.hdr.magic = 'ni1';
0363     elseif any(strcmp(suf,{'.nii','.NII'}))
0364         val.offset = max(sval.offset,352);
0365         obj.hdr.magic = 'n+1';
0366     else
0367         error(['Unknown filename extension (' suf ').']);
0368     end;
0369     val.offset        = (ceil(val.offset/16))*16;
0370     obj.hdr.vox_offset = val.offset;
0371 
0372     obj.hdr.dim(2:(numel(sz)+1)) = sz;
0373     nd = max(find(obj.hdr.dim(2:end)>1));
0374     if isempty(nd), nd = 3; end;
0375     obj.hdr.dim(1)   = nd;
0376     obj.hdr.datatype = sval.dtype;
0377     obj.hdr.bitpix   = d.size*8;
0378     if ~isempty(sval.scl_slope), obj.hdr.scl_slope = sval.scl_slope; end;
0379     if ~isempty(sval.scl_inter), obj.hdr.scl_inter = sval.scl_inter; end;
0380     obj.dat          = val;
0381 else
0382     error('"raw" must be of class "create_file_array"');
0383 end;
0384 return;
0385 
0386 function ok = valid_fields(val,allowed)
0387 if isempty(val), ok = false; return; end;
0388 if ~isstruct(val),
0389     error(['Expecting a structure, not a ' class(val) '.']);
0390 end;
0391 fn = fieldnames(val);
0392 for ii=1:length(fn),
0393     if ~any(strcmpi(fn{ii},allowed)),
0394         fprintf('Allowed fieldnames are:\n');
0395         for i=1:length(allowed), fprintf('    %s\n', allowed{i}); end;
0396         error(['"' fn{ii} '" is not a valid fieldname.']);
0397     end
0398 end
0399 ok = true;
0400 return;