Commit 3262bfb8 authored by Robert Henzel's avatar Robert Henzel

Merge branch 'master' into bassystExport

parents fd6be789 82cd51b1
...@@ -7,6 +7,7 @@ HTML ...@@ -7,6 +7,7 @@ HTML
Thumbs.db Thumbs.db
.DS_Store .DS_Store
*.asv *.asv
*.m~
*.wav *.wav
*.xlsx *.xlsx
*.docx *.docx
......
...@@ -788,15 +788,22 @@ classdef itaHRTF < itaAudio ...@@ -788,15 +788,22 @@ classdef itaHRTF < itaAudio
end end
function display(this) function display(this)
this.displayLineStart if numel(this) == 0
this.disp disp('****** nothing to do, empty object ******')
elseif numel(this) > 1
dir = num2str(this.nDirections,5); disp(['size(' inputname(1) ') = [' num2str(size(this)) ']; (for full display, pick a single instance)']);
stringD = [dir ' Directions (Type = ' this.mTF_type ')']; else
this.displayLineStart
this.disp
dir = num2str(this.nDirections,5);
stringD = [dir ' Directions (Type = ' this.mTF_type ')'];
middleLine = this.LINE_MIDDLE;
middleLine(3:(2+length(stringD))) = stringD;
fprintf([middleLine '\n']);
end
middleLine = this.LINE_MIDDLE;
middleLine(3:(2+length(stringD))) = stringD;
fprintf([middleLine '\n']);
end end
function disp(this) function disp(this)
......
...@@ -6,7 +6,7 @@ function [ data, samplerate, metadata ] = dfitaHRIRDAFFDataFunc( alpha, beta, it ...@@ -6,7 +6,7 @@ function [ data, samplerate, metadata ] = dfitaHRIRDAFFDataFunc( alpha, beta, it
% DAFF requires data alignment by multiple of 4 % DAFF requires data alignment by multiple of 4
nResidual = mod( hrtf.nSamples, 4 ); nResidual = mod( hrtf.nSamples, 4 );
data = [ hrtf.timeData', zeros( hrtf.nChannels, 4 - nResidual ) ]; data = [ hrtf.timeData', zeros( hrtf.nChannels, mod(4 - nResidual,4) ) ];
metadata = []; metadata = [];
end end
...@@ -44,26 +44,33 @@ end ...@@ -44,26 +44,33 @@ end
theta_start_deg = rad2deg( min( this.channelCoordinates.theta ) ); theta_start_deg = rad2deg( min( this.channelCoordinates.theta ) );
theta_end_deg = rad2deg( max( this.channelCoordinates.theta ) ); theta_end_deg = rad2deg( max( this.channelCoordinates.theta ) );
theta_num_elements = size( unique( this.channelCoordinates.theta ), 1 ); theta_num_elements = size( uniquetol( this.channelCoordinates.theta ), 1 );
phi_start_deg = rad2deg( min( mod( this.channelCoordinates.phi, 2*pi ) ) ); phi_start_deg = rad2deg( min( mod( this.channelCoordinates.phi, 2 * pi ) ) );
phi_end_deg = rad2deg( max( mod( this.channelCoordinates.phi, 2*pi ) ) ); phi_end_deg = rad2deg( max( mod( this.channelCoordinates.phi, 2 * pi ) ) );
phi_num_elements = size( unique( this.channelCoordinates.phi ), 1 ); phi_num_elements = size( uniquetol( this.channelCoordinates.phi ), 1 );
assert( phi_num_elements ~= 0 ); assert( phi_num_elements ~= 0 );
alphares = ( phi_end_deg - phi_start_deg ) / phi_num_elements; % phi end does not cover entire circle in this case alphares = ( phi_end_deg - phi_start_deg ) / phi_num_elements; % phi end does not cover entire circle in this case
alphares_full_circle = ( phi_end_deg - phi_start_deg ) / ( phi_num_elements - 1 ); % phi end does not cover entire circle in this case alphares_full_circle = ( phi_end_deg - phi_start_deg ) / ( phi_num_elements - 1 ); % phi end does not cover entire circle in this case
if phi_end_deg + alphares_full_circle >= 360.0 if phi_end_deg + alphares_full_circle >= 360.0
alpharange = [ phi_start_deg ( phi_end_deg + alphares_full_circle ) ]; % Account for full circle alpharange = [ phi_start_deg 360 ]; % Account for full circle and force end of range to 360 deg
alphares = alphares_full_circle; alphares = alphares_full_circle;
else else
alpharange = [ phi_start_deg phi_end_deg ]; alpharange = [ phi_start_deg phi_end_deg ];
end end
assert( alpharange( 1 ) >= 0.0 )
assert( alpharange( 2 ) <= 360.0 )
assert( theta_num_elements ~= 0 ); assert( theta_num_elements ~= 0 );
betares = ( theta_end_deg - theta_start_deg ) / ( theta_num_elements - 1 ); % phi end does not cover entire circle betares = ( theta_end_deg - theta_start_deg ) / ( theta_num_elements - 1 ); % phi end does not cover entire circle
betarange = 180 - [ theta_start_deg theta_end_deg ]; % Flip poles (DAFF starts at south pole) betarange = 180 - [ theta_start_deg theta_end_deg ]; % Flip poles (DAFF starts at south pole)
assert( betarange( 2 ) >= 0.0 )
assert( betarange( 1 ) <= 180.0 )
%% Assemble metadata %% Assemble metadata
metadata = daffv17_add_metadata( metadata, 'Generation script', 'String', 'writeDAFFFile.m' ); metadata = daffv17_add_metadata( metadata, 'Generation script', 'String', 'writeDAFFFile.m' );
......
...@@ -41,7 +41,7 @@ HRTF_find.play_gui(pinkNoise); ...@@ -41,7 +41,7 @@ HRTF_find.play_gui(pinkNoise);
%% Binaural parameters %% Binaural parameters
ITD = slicePhi.ITD; % different methods are available: see method in itaHRTF ITD = slicePhi.ITD; % different methods are available: see method in itaHRTF
ILD = slicePhi.ILD; %ILD = slicePhi.ILD;
%% Modifications %% Modifications
% calculate DTF % calculate DTF
...@@ -65,8 +65,12 @@ HRTF_interp = HRTF_sphere.interp(coordI); ...@@ -65,8 +65,12 @@ HRTF_interp = HRTF_sphere.interp(coordI);
nameDaff_file = 'HRTF_sphere.daff'; nameDaff_file = 'HRTF_sphere.daff';
HRTF_sphere.writeDAFFFile(nameDaff_file); HRTF_sphere.writeDAFFFile(nameDaff_file);
HRTF_daff = itaHRTF('daff',nameDaff_file); %HRTF_daff = itaHRTF('daff',nameDaff_file);
nameDaff_file2 = 'C:\Users\bomhardt\Documents\ITA-Toolbox\applications\TODOunfinished\test_marcia\DATA\HRTF_ITAKopf_Nov2013_MartinPollow\ITA-Kunstkopf_HRIR_Mess01_D180_1x1_256.daff'; nameDaff_file2 = 'yourHRTF.daff';
HRTF_daff2 = itaHRTF('daff',nameDaff_file2); if ~strcmp(nameDaff_file2,'yourHRTF.daff')
HRTF_daff2.plot_freqSlice HRTF_daff2 = itaHRTF('daff',nameDaff_file2);
\ No newline at end of file HRTF_daff2.plot_freqSlice
else
ita_disp('use an existing daff-file')
end
\ No newline at end of file
...@@ -14,6 +14,7 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -14,6 +14,7 @@ classdef itaHpTF_Audio < itaHpTF
% Audio Engineering Society Convention 130, May 2011) % Audio Engineering Society Convention 130, May 2011)
% normalized (Normalization of HpTF) % normalized (Normalization of HpTF)
% smoothing (Smoothing bandwidth has to be defined in fractions of octaves, see also ita_smooth) % smoothing (Smoothing bandwidth has to be defined in fractions of octaves, see also ita_smooth)
% mic (measured microphone transfer funcions - not inverted: length does not matter; will be adapted)
% %
% itaHpTF_Audio Methods (private): % itaHpTF_Audio Methods (private):
% HP_equalization If this.mic is not set, it will be unused % HP_equalization If this.mic is not set, it will be unused
...@@ -30,7 +31,9 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -30,7 +31,9 @@ classdef itaHpTF_Audio < itaHpTF
m_fUpper = 18000; m_fUpper = 18000;
mMethod = 'mSTD'; mMethod = 'mSTD';
mNormalized = true; mNormalized = true;
mGainComp = true;
mSmoothing = 1/6; mSmoothing = 1/6;
mSelectMeas = 1:8; % select the channels which should be used for the HpTF
end end
properties(Dependent = true, Hidden = false) properties(Dependent = true, Hidden = false)
...@@ -38,8 +41,10 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -38,8 +41,10 @@ classdef itaHpTF_Audio < itaHpTF
fLower = 100; fLower = 100;
fUpper = 18000; fUpper = 18000;
method = 'mSTD'; method = 'mSTD';
normalized = true; normalized = true; % normalization of the signal
gainComp = true; % microphone compensation
smoothing = 1/6; smoothing = 1/6;
selectMeas = 1:8; % select specific measurements
end end
properties (Dependent = true, SetAccess = private) properties (Dependent = true, SetAccess = private)
...@@ -122,6 +127,13 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -122,6 +127,13 @@ classdef itaHpTF_Audio < itaHpTF
function normalized = get.normalized(this) function normalized = get.normalized(this)
normalized = this.mNormalized; end normalized = this.mNormalized; end
function comp = get.gainComp(this)
comp = this.mGainComp; end
function sel = get.selectMeas(this)
sel = this.mSelectMeas; end
function smoothing = get.smoothing(this) function smoothing = get.smoothing(this)
smoothing = this.mSmoothing; end smoothing = this.mSmoothing; end
%% ................................................................ %% ................................................................
...@@ -151,6 +163,16 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -151,6 +163,16 @@ classdef itaHpTF_Audio < itaHpTF
this = HP_equalization(this); this = HP_equalization(this);
end end
function this = set.selectMeas(this,ch)
this.mSelectMeas = ch;
this = HP_equalization(this);
end
function this = set.gainComp(this,comp)
this.mGainComp = comp;
this = HP_equalization(this);
end
function this = set.smoothing(this,smoothing) function this = set.smoothing(this,smoothing)
this.mSmoothing = smoothing; this.mSmoothing = smoothing;
this = HP_equalization(this); this = HP_equalization(this);
...@@ -162,11 +184,12 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -162,11 +184,12 @@ classdef itaHpTF_Audio < itaHpTF
% SET PRIVATE % SET PRIVATE
%.................................................................. %..................................................................
function this = set.init(this,var) function this = set.init(this,var)
this.nameHP = var.nameHP; this.nameHP = var.nameHP;
this.nameMic = var.nameMic; this.nameMic = var.nameMic;
this.nameSubj = var.nameSubj; this.nameSubj = var.nameSubj;
this.repeat = var.repeat; this.repeat = var.repeat;
this.mTF.time = var.time; this.mTF.time = var.time;
this.selectMeas = 1:var.repeat;
this = HP_equalization(this); this = HP_equalization(this);
end end
...@@ -175,11 +198,8 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -175,11 +198,8 @@ classdef itaHpTF_Audio < itaHpTF
% Audio Engineering Society Convention 130, May 2011 % Audio Engineering Society Convention 130, May 2011
tWin = this.TF.trackLength; % crop HPTF tWin = this.TF.trackLength; % crop HPTF
if this.mic.dimensions == 2 % mic min phase + extension measSel = sort([2*this.selectMeas-1 2*this.selectMeas]);
minMic = ita_minimumphase(this.mic); TF = this.TF.ch(measSel);
[mic, TF] = ita_extend_dat(minMic,this.TF);
else TF = this.TF;
end
%init %init
thisEQ = this; thisEQ = this;
...@@ -200,14 +220,19 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -200,14 +220,19 @@ classdef itaHpTF_Audio < itaHpTF
otherwise otherwise
error('Unknown type'); error('Unknown type');
end end
if this.mic.dimensions == 2
Rec = ita_multiply_spk(Rec,mic); if this.mic.dimensions == 2 % Compensation of mic
minMic = ita_minimumphase(this.mic);
RecM = ita_convolve(Rec,minMic.ch(cdx));
RecM.fftDegree = TF.fftDegree;
else
RecM = Rec;
end end
%% Short Filter with no correction for low freqs and minimun phase %% Short Filter with no correction for low freqs and minimum phase
aux = max(abs(Rec.freqData),[],2); aux = max(abs(RecM.freqData),[],2);
% find first maximum and truncate low freq correction at this point % find first maximum and truncate low freq correction at this point
idxDF = Rec.freq2index([this.fLower this.fLower*1.5 ]); idxDF = RecM.freq2index([this.fLower this.fLower*1.5 ]);
d_aux = diff(aux(idxDF(1):idxDF(2))); d_aux = diff(aux(idxDF(1):idxDF(2)));
idx = find(diff(sign(d_aux)) ~= 0,1,'first'); % Bruno style... idx = find(diff(sign(d_aux)) ~= 0,1,'first'); % Bruno style...
aux(1:idxDF(1)+idx+1) = aux(idxDF(1)+idx+2); aux(1:idxDF(1)+idx+1) = aux(idxDF(1)+idx+2);
...@@ -227,10 +252,18 @@ classdef itaHpTF_Audio < itaHpTF ...@@ -227,10 +252,18 @@ classdef itaHpTF_Audio < itaHpTF
this_min = ita_minimumphase(ita_time_shift(HpTF,tWin/2)); this_min = ita_minimumphase(ita_time_shift(HpTF,tWin/2));
this_win = ita_time_window(this_min,[tWin*0.99,tWin],'time','crop'); this_win = ita_time_window(this_min,[tWin*0.99,tWin],'time','crop');
if this.normalized, if this.gainComp % compensate gains from mics (equalize)
this_norm = ita_normalize_dat(this_win); idxDF = this_win.freq2index(this.fLower);
chGain = 20*log10(abs(this_win.freqData(idxDF,:)));
this_comp = this_win;
this_comp.freqData = bsxfun(@times,this_win.freqData,10.^(-chGain/20));
else, this_comp = this_win;
end
if this.normalized
this_norm = ita_normalize_dat(this_comp);
thisEQ.timeData = this_norm.timeData; thisEQ.timeData = this_norm.timeData;
else thisEQ.timeData = this_win.timeData; else, thisEQ.timeData = this_comp.timeData;
end end
if ~isempty(this.savePath) if ~isempty(this.savePath)
......
...@@ -88,6 +88,7 @@ classdef itaHpTF_MS < itaHpTF ...@@ -88,6 +88,7 @@ classdef itaHpTF_MS < itaHpTF
% Measure Right side % Measure Right side
MS.inputChannels = chIn(2); MS.inputChannels = chIn(2);
MS.outputChannels = chOut(2); MS.outputChannels = chOut(2);
commandwindow
resultR = MS.run; resultR = MS.run;
resultR.channelNames{1} = strR; resultR.channelNames{1} = strR;
...@@ -108,7 +109,7 @@ classdef itaHpTF_MS < itaHpTF ...@@ -108,7 +109,7 @@ classdef itaHpTF_MS < itaHpTF
disp(['Measurement ',num2str(idxM),'/',num2str(this.repeat),' finished.']) disp(['Measurement ',num2str(idxM),'/',num2str(this.repeat),' finished.'])
disp('Put headphones off and on again. Press any key to continue.') disp('Put headphones off and on again. Press any key to continue.')
pause pause
pause(2) pause(1)
else else
commandwindow commandwindow
fprintf('\nMEASUREMENT DONE!\n') fprintf('\nMEASUREMENT DONE!\n')
...@@ -118,6 +119,7 @@ classdef itaHpTF_MS < itaHpTF ...@@ -118,6 +119,7 @@ classdef itaHpTF_MS < itaHpTF
MS.inputChannels = chIn; MS.inputChannels = chIn;
HpTF = itaHpTF_Audio(this); HpTF = itaHpTF_Audio(this);
HpTF.selectMeas = 1:this.repeat;
assignin('base', ['HpTF_' this.nameSubj], HpTF) assignin('base', ['HpTF_' this.nameSubj], HpTF)
end end
......
classdef itaSerialDeviceInterface < handle
% Wrapper around a serial object. Contains initialization code and wraps functions. Is a
% Singleton so that outside objects can access the serial object. If communication is to
% take place over another interface (e.g. CAN), reuse interface of this class by abstracting
% away the specific stuff in here, then introducing an abstract base class itaMotorInterface,
% then inherit from it.
properties
mSerialObj;
comPort = 'uninitialized';
portOpen = false;
baudRate = 19200;
databits = 8;
stopbits = 1;
OutputBufferSize = 3072;
Terminator = 13;
BytesAvailableFcnMode = 'Terminator';
bytesAvailable;
end
% Constructor, Destructor ------------------------------------------------------------------
methods %(Access = private) TODO : make private again to prevent ordinary instantiation?
function this = itaSerialDeviceInterface(varargin)
% config
this.comPort = ita_preferences('movtecComPort');
if strcmpi(this.comPort,'noDevice')
ita_verbose_info('Please select a COM-Port in ita_preferences (Using the Movtec-Port!)',0); % TODO : messaging
ita_preferences;
return;
end
% init serial object
this.initComPort();
end
function delete(this)
% TODO : check if serial object is open first?
fclose( this.mSerialObj);
end
end
% Public interface ------------------------------------------------------------------------
methods
function setBaudRate(this,br)
if ~any(this.validBaudRates == br)
ita_verbose_info(sprintf('Invalid baudRate selected: %d, Valid baud rates: %s',this.baudRate, num2str(this.validBaudRates)),0); % TODO : messaging
return;
end
this.mSerialObj.baudRate = br;
end
function ret = query(this, msg)
if(this.portOpen == false)
this.portOpen = true;
ret = query(this.mSerialObj, sprintf(msg));
this.portOpen = false;
else
ret = query(this.mSerialObj, sprintf(msg));
end
end
function sendAsynch(this,msg)
% disp(sprintf('send: %s',msg));
fwrite(this.mSerialObj, sprintf(msg));
end
function ret = recvAsynch(this)
ret = fgetl(this.mSerialObj);
% disp(sprintf('rec: %s',ret));
end
function avail = BytesAvailable(this)
avail = this.mSerialObj.bytesAvailable;
end
end
% Singleton ------------------------------------------------------------------------------------
methods (Static)
% Singleton access function
function instance = getInstance
persistent localInstance
if isempty(localInstance) || ~isvalid(localInstance)
localInstance = itaSerialDeviceInterface;
end
instance = localInstance;
end
end
% Getters n Setters ----------------------------------------------------------------------------
methods
function set.baudRate(this, value)
this.mSerialObj.BaudRate = value;
this.baudRate = value;
end
function set.portOpen(this,value)
if(this.portOpen ~= true)
% TODO : if this is the first time that the serial interface
% is accessed, check if com port is correct?
fopen(this.mSerialObj);
this.portOpen = true;
end
end
function set.comPort(this,value)
this.comPort = value;
this.initComPort();
end
function val = get.bytesAvailable(this)
val = this.mSerialObj.BytesAvailable;
end
end
% Private functions ----------------------------------------------------------------------------
methods (Access = private)
function initComPort(this)
insts = instrfind; %show existing terminals using serial interface
if ~isempty(insts)
aux = strfind(insts.Name,this.comPort);
if numel(aux) == 1
aux = {aux};
end
for idx = 1:numel(aux)
if ~isempty(aux{idx})
delete(insts(idx)); %delete used serial ports
end
end
end
this.mSerialObj = serial(this.comPort,'Baudrate',this.baudRate,'Databits',this.databits,'Stopbits',this.stopbits,'OutputBufferSize',this.OutputBufferSize);
this.mSerialObj.Terminator = this.Terminator;
this.mSerialObj.BytesAvailableFcnMode = this.BytesAvailableFcnMode;
end
end
end
\ No newline at end of file
classdef itaMotorControl < itaHandle
%ITAMOTORCONTROL Summary of this class goes here
% Detailed explanation goes here
properties(Access = protected, Hidden = true)
mSerialObj; % the serial connection
end
properties
comPort = ita_preferences('movtecComPort');
baudrate = 9600;
databits = 8;
stopbits = 1;
OutputBufferSize = 3072;
wait = true; % Status - do we wait for motors to reach final position? -> Set by prepare-functions!
end
methods(Abstract)
% basic motor functions
this = init(this);
stopAllMotors(this);
setWait(this,value);
% isInitialized(this);
% basic moves: requires execution to halt while something is moving
this = reference(this);
this = moveTo(this,targetPosition);
ret = prepareForContinuousMeasurement(this);
startContinuousMoveNow(this);
end
methods(Abstract, Hidden=true)
displayMotors(this)
end
end
classdef itaMotorNanotec < itaHandle
properties(Access = protected, Hidden = true)
motorID = [];
motorName = [];
% mInUse = false;
old_position = itaCoordinates(1);
mIsReferenced = false;
mIsInit = false;
motorLimits = [];
end
properties
end
methods(Abstract)
% basic motor functions
this = init(this);
stop(this);
active = isActive(this);
setActive(this,value);
% isInitialized(this);
getStatus(this);
getMotorID(this);
getMotorName(this);
% basic moves: requires execution to halt while something is moving
this = moveToReferencePosition(this);
this = startMoveToPosition(this);
ret = prepareMove(this,position,varargin);
end
end
\ No newline at end of file
...@@ -553,7 +553,7 @@ classdef itaOptitrack < handle ...@@ -553,7 +553,7 @@ classdef itaOptitrack < handle
if Optitrack_obj.singleShot if Optitrack_obj.singleShot
Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once
else else
Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,'ExecutionMode','fixedRate','BusyMode','drop'); Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,'ExecutionMode','fixedSpacing','BusyMode','drop');
end end
else else
Optitrack_obj.correctRowIdx = 1; % reset row counter Optitrack_obj.correctRowIdx = 1; % reset row counter
...@@ -561,7 +561,7 @@ classdef itaOptitrack < handle ...@@ -561,7 +561,7 @@ classdef itaOptitrack < handle
Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'ExecutionMode','singleShot'); % execute timer callback once
else else
Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,... Optitrack_obj.timerData = timer('TimerFcn',{@Optitrack_obj.TimerCallback},'Period',round(1/16/Optitrack_obj.frameRate*1000)/1000,...
'ExecutionMode','fixedRate','BusyMode','drop','TasksToExecute',32*Optitrack_obj.numFrames); % execute timer 16*Optitrack_obj.numFrames times but abort as soon as Optitrack_obj.rigidBodyLogData.data is filled 'ExecutionMode','fixedSpacing','BusyMode','drop','TasksToExecute',32*Optitrack_obj.numFrames); % execute timer 16*Optitrack_obj.numFrames times but abort as soon as Optitrack_obj.rigidBodyLogData.data is filled
end end
end end
......
...@@ -52,8 +52,15 @@ if abs(dot(v,u,2)) > 1e-5 ...@@ -52,8 +52,15 @@ if abs(dot(v,u,2)) > 1e-5
end end
% normalize view/up vectors % normalize view/up vectors
v = normr(v); [~,colv]=size(v);
u = normr(u); [~,colu]=size(u);
if (colv == 1)
v(~isnan(v(:,1)),:) = v(~isnan(v(:,1)),:) ./ abs(v(~isnan(v(:,1)),:));
u(~isnan(u(:,1)),:) = u(~isnan(u(:,1)),:) ./ abs(u(~isnan(u(:,1)),:));
else
v = sqrt( ones ./ (sum((v(~isnan(v(:,1)),:).*v(~isnan(v(:,1)),:))')) )' * ones(1,colv).*v(~isnan(v(:,1)),:);
u = sqrt( ones ./ (sum((u(~isnan(u(:,1)),:).*u(~isnan(u(:,1)),:))')) )' * ones(1,colu).*u(~isnan(u(:,1)),:);
end
% calculate side vector % calculate side vector
s = cross(v, u); s = cross(v, u);
...@@ -65,7 +72,7 @@ qx = qw; ...@@ -65,7 +72,7 @@ qx = qw;
qy = qw; qy = qw;
qz = qw; qz = qw;
for idx = 1:vec_ent; for idx = 1:vec_ent