a = csvread('original pixel data.csv');
a = a*0.045823; % convert voxels to mm
origin = a(1,:);% use point 1 for origin
num = size(a);
for i=1:num(1)
a(i,:) = a(i,:) - origin; %remove offset from point 1
end

% establish origin and points used to define axes with offset data
origin = a(1,:);% use point 1
zpnt = a(31,:); % use point 31
xpnt = a(6,:);  % use point 6

zaxis = origin - zpnt; % vector from pnt 31 to pnt 1
xaxis = xpnt - origin; % vector from pnt 31 to pnt 1

znorm = zaxis/norm(zaxis); %Normalize z vector
xtmpnorm = xaxis/norm(xaxis); %Normalize xtemp vector

yaxis = cross(znorm,xtmpnorm); % find y axis
ynorm = yaxis/norm(yaxis);%Normalize y vector

xaxis = cross(ynorm,znorm); % find x axis. since x and z may not have been perpindicular
xnorm = xaxis/norm(xaxis);%Normalize x vector.  Just to be sure

R = [xnorm' ynorm' znorm']; % create rotation matrix
calpoints = (R'*a')'; % transform coordinates to localized CS

%gather information for file header
Name = input('What is the name of the person who created with file? ', 's');
device = input('What is the name of the device that was being calibrated? ', 's');

% create header file
CalFileHeader = {'Date/Time', datestr(now);
                 'Created By', Name;
                 'Calibration Device', device;
                 'Number of Points', num(1)};

% create units
               units = 'mm';
                
% write calibration file
save('DevACalFile_full.mat', 'CalFileHeader', 'calpoints', 'units');

