vb_read_curry_surface

 Read surface data from Curry format file
  [F,V,xx] = vb_read_curry_surface(curryfile, analyzefile)

 --- Input parameters
 curryfile
 analyzefile
 --- Output
 V  : vertex coordinate   : NV x 3 
 F  : surface patch index : NF x 3
 xx : normal vector       : NV x 3 

 Curryで作成されたサーフェイス情報(*.s00-s99)をロードする。
 ただし、Curryの起動前にstartup.cfgを変更し、SAVE_PACKEDを
 0に設定しておくこと。
 --- DEBUG info
  [F,V,xx,Ninfo] = vb_read_curry_surface(curryfile, analyzefile)
 Ninfo : Debug info

 2004-01-30 Taku Yoshioka
 2004-02-12 Taku Yoshioka
 - 出力引数の順序を変更(Fが先になるようにした)
 2005-07-06 Taku Yoshioka
 - 右手系SPM-ANALYZE座標系で出力する
 2006-10-26 M.Sato
  ANALYZEファイルの引数を削除
  vb_create_curry_patch を高速化
 2007-06-25 M.Sato
  任意の大きさの画像に対応するためANALYZEファイルの引数を復活
 2009-12-15 M.Sato
  voxcel サイズが1mm以外の画像に対応
---

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [F,V,xx,Ninfo] = vb_read_curry_surface(curryfile, analyzefile)
0002 % Read surface data from Curry format file
0003 %  [F,V,xx] = vb_read_curry_surface(curryfile, analyzefile)
0004 %
0005 % --- Input parameters
0006 % curryfile
0007 % analyzefile
0008 % --- Output
0009 % V  : vertex coordinate   : NV x 3
0010 % F  : surface patch index : NF x 3
0011 % xx : normal vector       : NV x 3
0012 %
0013 % Curryで作成されたサーフェイス情報(*.s00-s99)をロードする。
0014 % ただし、Curryの起動前にstartup.cfgを変更し、SAVE_PACKEDを
0015 % 0に設定しておくこと。
0016 % --- DEBUG info
0017 %  [F,V,xx,Ninfo] = vb_read_curry_surface(curryfile, analyzefile)
0018 % Ninfo : Debug info
0019 %
0020 % 2004-01-30 Taku Yoshioka
0021 % 2004-02-12 Taku Yoshioka
0022 % - 出力引数の順序を変更(Fが先になるようにした)
0023 % 2005-07-06 Taku Yoshioka
0024 % - 右手系SPM-ANALYZE座標系で出力する
0025 % 2006-10-26 M.Sato
0026 %  ANALYZEファイルの引数を削除
0027 %  vb_create_curry_patch を高速化
0028 % 2007-06-25 M.Sato
0029 %  任意の大きさの画像に対応するためANALYZEファイルの引数を復活
0030 % 2009-12-15 M.Sato
0031 %  voxcel サイズが1mm以外の画像に対応
0032 %---
0033 %
0034 % Copyright (C) 2011, ATR All Rights Reserved.
0035 % License : New BSD License(see VBMEG_LICENSE.txt)
0036 
0037 fprintf('--- Read curry surface data \n');
0038 fprintf('File = %s\n', curryfile);
0039 
0040 if ~exist(curryfile,'file'), error('No curry file');end;
0041 
0042 fp = fopen(curryfile,'r');
0043 
0044 if fp < 0, error('Can not open file'), end;
0045 
0046 %
0047 % 頂点数('LOCATION_LIST START'の後ろの情報)
0048 %
0049 while 1
0050   str = fgets(fp);
0051   if ~isempty(findstr(str,'LOCATION_LIST START'))
0052     while 1
0053       str = fgets(fp);
0054       if ~isempty(findstr(str,'LIST_NR_ROWS'))
0055         ix = findstr(str,'=');
0056         NV0 = str2num(str(ix+1:length(str)));
0057         fprintf('Number of faces: %d\n',NV0/3);
0058         break;
0059       end
0060 
0061       if feof(fp)
0062         error(['Entry ''LOCATION_LIST START'' is not closed.']);
0063       end
0064     end
0065   end
0066   
0067   if exist('NV0')
0068     break;
0069   end
0070   
0071   if feof(fp)
0072     error(['There is no entry ''LOCATION_LIST START'' ' ...
0073        'in file ' curryfile '.']);
0074   end
0075 end
0076 
0077 disp('--- Read vertex position');
0078 
0079 %
0080 % 頂点座標開始点('LOCATION_LIST START_LIST'まで読み飛ばす)
0081 %
0082 while 1
0083   str = fgets(fp);
0084   if ~isempty(findstr(str,'LOCATION_LIST START_LIST'))
0085     break;
0086   end
0087 
0088   if feof(fp)
0089     error(['There is no entry ''LOCATION_LIST START_LIST'' ' ...
0090        'in file ' curryfile '.']);
0091   end
0092 end
0093 
0094 %
0095 % --- Read Curry surface
0096 %
0097 FV = (fscanf(fp,'%u %u %u',[3 NV0]))';
0098 fclose(fp);
0099 
0100 %
0101 % --- Create patch index from Curry vertex coordinate
0102 %
0103 tic;
0104 [V,F,ix,Ninfo] = vb_create_curry_patch(FV);
0105 vb_ptime(toc)
0106 
0107 % ---  Debug
0108 %fprintf('Npos=%d\n',size(FV,1))
0109 
0110 if Ninfo.NFcnt ~= size(FV,1)
0111     fprintf('# of triangle vertex = %d\n',size(FV,1))
0112     fprintf('# of processed vertex\n',Ninfo.NFcnt)
0113     error('# of processed vertex is different')
0114 end
0115 if Ninfo.NVcnt ~= size(V,1)
0116     fprintf('# of vertex = %d\n',size(V,1))
0117     fprintf('# of processed vertex\n',Ninfo.NVcnt)
0118     error('# of processed vertex is different')
0119 end
0120 
0121 fprintf('# of vertex = %d\n',size(V,1))
0122 fprintf('# of patch  = %d\n',size(F,1))
0123 
0124 %
0125 % 法線情報(NORMALS_LIST)
0126 %
0127 disp('--- Read normal vector');
0128 fp = fopen(curryfile,'r');
0129 
0130 while 1
0131   str = fgets(fp);
0132   if ~isempty(findstr(str,'NORMALS_LIST START_LIST'))
0133     break;
0134   end
0135 
0136   if feof(fp)
0137     error(['There is no entry ''NORMALS_LIST START_LIST'' ' ...
0138        'in file ' curryfile '.']);
0139   end
0140 end
0141 
0142 xx = (fscanf(fp,'%d %d %d',[3 NV0]))';
0143 fclose(fp);
0144 
0145 % outward normal direction
0146 xx = xx(ix,:);
0147 xx = -1*xx./repmat(sqrt(sum(xx.^2,2)),[1 3]);
0148 
0149 fprintf('# of normal = %d\n',size(xx,1))
0150 
0151 %%%%% ---- DEBUG
0152 %return
0153 %%%%% ---- DEBUG
0154 
0155 %
0156 % --- Coordinate transformation from Curry to SPM-right
0157 %
0158 
0159 % --- Curry mm coordinate
0160 % [Curry]
0161 % Xc:Right(1)->Left(256)
0162 % Yc:Front(1)->Back(256)
0163 % Zc:Bottom(1)->Top(256)
0164 % --- SPM mm coordinate
0165 % Xs: Left(-NX/2)   -> Right(NX/2)
0166 % Ys: Back(-NY/2)   -> Front(NY/2)
0167 % Zs: Bottom(-NZ/2) -> Top(NZ/2)
0168 % --- adjustment point : XY-centor & Z-Top
0169 % RL-mid : Xc = 128 : Xs = 0
0170 % FB-mid : Yc = 128 : Ys = 0
0171 % Top    : Zc = 256 : Zs = NZ/2
0172 % ---
0173 % Xs = - (Xc - NX/2)
0174 % Ys = - (Yc - NY/2)
0175 % Zs = Zc + NZ/2 - 256
0176 
0177 [Vdim, Vsize] = analyze_hdr_read(analyzefile);
0178 
0179 X0 = 128;
0180 Y0 = 128;
0181 Z0 = 256;
0182 
0183 % adjust Curry origin
0184 V(:,1) = V(:,1) - X0;
0185 V(:,2) = V(:,2) - Y0;
0186 V(:,3) = V(:,3) - Z0;
0187 
0188 % scale transformation to mm
0189 scale  = Vsize(2);
0190 %scale  = min(Vsize);
0191 
0192 V  = scale*V;
0193 
0194 % adjust SPM origin
0195 X0 = 0;
0196 Y0 = 0;
0197 Z0 = Vdim(3)*Vsize(3)/2;
0198 
0199 % Curry mm to SPM mm coordinate
0200 V(:,1) = X0 - V(:,1);
0201 V(:,2) = Y0 - V(:,2);
0202 V(:,3) = Z0 + V(:,3);
0203 
0204 xx(:,1) = -xx(:,1);
0205 xx(:,2) = -xx(:,2);
0206 xx(:,3) =  xx(:,3);
0207 
0208 V  = V./1000;    % mm -> m
0209 
0210 solid_angle  = vb_solid_angle_check(V,F);
0211 
0212 fprintf('Closed surface index = %f ( = 1)\n',solid_angle);