Octave functions for generating Bezier yacht hull curves using Letcher's method

BezierHull.m

% Generate yacht hull lines using Letcher's method with all curves defined
% analytically using Bezier curves
%                                      
%                                      
%                                      
% Input(s)                             
%   fName: Name of Excel file containing Bezier points                                  
%                                      
% Output(s)                            
%   BezStruct: 
%    BezCurves: Structure of Bezier curves for each of the 6 fundamental curves
%               Sheer, Freeboard, Profile, A (homotopy), eta1 and eta2
%    sections:  Section curves derived from the Bezier functions at nSections 
%               intervals along the hull
%    functions: Symbolically defined functions for each section (eta), and a 
%               final set for the transom, (trans_x, trans_y, trans_z)
%    Boat parameters:
%          LOA: Length overall
%         Beam: Maximum beam
%        Draft: Maximum draft
%      Midship: Distance from bow of maximum beam
%
% <fName>.txt used as input to DelftShip is written to the local directory
%                
%                                                 
%                                      
% Example:                             
% pkg load io
% pkg load symbolic 
% BezStruct = BezierHull('Letcher curves.xlsx');   
% Start DELFTship, click on "New Project" (Ctrl-N)
% Select a blank project (white rectangle in toolstrip) and "Accept"
% Click on "Open" down arrow -> "Import" -> "Surface"                                
%                                      
% See also:                            
%                                      
%                                      
% Dependencies: Packages: io, symbolic
%               Functions: scale, symBezFunc, funcSolve, write2DelftShip
                                       
%                                      
%                                      
% Written by: John Peach 26-Sep-2020   
% Marakesh Sailing Analysis and Design 
%                                      
% Revisions:                           
                                       
                                       
function BezStruct = BezierHull(fName)
    
    % Distances between sections and waterlines
    dx = 1;     % Section distance
    dz = 0.25;  % Waterline distance
  
    % Read spreadsheet data
    disp('Reading Excel file')
    [~,~,df] = xlsread(fName);
    
    % Separate into curve data
    crvNames = {'Sheer','Freeboard','Profile','A','eta1','eta2'};
    BezCurves = repmat(struct('Curve',[],'x',[],'y',[]),6,1);
    crvIdx = zeros(1,6);
    for k = 1:size(df,1)
        if ~isempty(df{k,1})
            i = cellfind(crvNames,df{k,1});
            if ~isempty(i)
                BezCurves(i).Curve = crvNames{i};
                crvIdx(i) = k;
            endif
        endif
    endfor
    crvIdx = [crvIdx size(df,1)+1];
    for k = 1:6
        BezCurves(k).x = [df{crvIdx(k):crvIdx(k+1)-1,2}];
        BezCurves(k).y = [df{crvIdx(k):crvIdx(k+1)-1,3}];
    endfor
        
    % Scale the A curve in the y dimension to [0,1]
    BezCurves(4).y = scale(BezCurves(4).y,[0, 1]);
    
    % Symbolic expressions for each function
    disp('Generating symbolic functions')
    p = symBezFunc(BezCurves,'Profile');  
    f = symBezFunc(BezCurves,'Freeboard');  
    s = symBezFunc(BezCurves,'Sheer');  
    A = symBezFunc(BezCurves,'A');  
    eta1 = symBezFunc(BezCurves,'eta1');  
    eta2 = symBezFunc(BezCurves,'eta2');
    
    % Positions of section curves. Assumes that the waterlines is at x = 0.
    x_max = p.x(1);
    x = 0:dx:x_max;       % No section curve drawn at waterlines
    nSections = numel(x);
    
    % Get values of t for each function at selected sections
    p.t = funcSolve(p.x,x);
    f.t = funcSolve(f.x,x);
    s.t = funcSolve(s.x,x);
    A.t = funcSolve(A.x,x);    
    
    % Waterline levels 
    z_min   = p.y(fminbnd(p.y,0,1));          % Minimum profile
    z_max   = f.y(fminbnd(@(t) -f.y(t),0,1)); % Maximum freeboard      
    zLevels = roundi(z_min,dz,'down'):dz:roundi(z_max,dz,'up');    
    nWaters = numel(zLevels); % Number of points along each section curve

    % Structures. Add one extra section for transom, two extra waterlines for
    % profile and sheer curves
    functions  = repmat(struct('eta',[]),nSections+1,1);
    sections   = repmat(struct('Pts',[]),nSections+1,1);
                      
    % Coordinates for points at each section  
    disp('Calculating section points')    
    for k = 1:nSections
        % Section sheer
        s.Section = s.y(s.t(k));
        % Section profile and freeboard
        p.Section = p.y(p.t(k));
        f.Section = f.y(f.t(k));
        % Section function at x(k) parameterized by t
        % eta1 defines the forward section curves, eta2 is for aft curves
        A.Section = A.y(A.t(k));  % Homotopy function
        % Section curves in y and z
        eta.y = @(t) s.Section * (A.Section * eta1.x(t) + ...
                    (1-A.Section) * eta2.x(t));
        eta.z = @(t) (f.Section-p.Section) * (A.Section * eta1.y(t) + ...
                     (1-A.Section) * eta2.y(t)) + p.Section;
        % Save section points
        t = [0 funcSolve(eta.z,zLevels) 1];
        i = find(~isnan(t));
        ni = numel(i);
        sections(k).Pts  = [x(k)*ones(1,ni); eta.y(t(i)); eta.z(t(i))]';
        % Section function
        functions(k).eta = eta;
        
     endfor 
     
     % Generate extra section for transom
     i = p.y(1) <= zLevels & zLevels <= f.y(1); % Transom waterlevels 
     trans_x = @(t) (f.x(1)-p.x(1))*t + p.x(1);
     trans_y = @(t) s.y(1) * eta2.x(t);
     trans_z = @(t) (f.y(1)-p.y(1))*eta2.y(t) + p.y(1);
     tt = funcSolve(trans_z,zLevels(i));
     transom = [trans_x(tt); trans_y(tt); trans_z(tt)]';  
     sections(nSections+1).Pts = transom;
     functions(end).eta.trans_x = trans_x;
     functions(end).eta.trans_y = trans_y;
     functions(end).eta.trans_z = trans_z;
     
     % Write data to DelftShip format
     disp('Writing points to DELFTship format file')
     write2DELFTship(sections,fName)
     
     % Display LOA, Beam, draft
     LOA = max([f.x(1) p.x(1)]);
     ts = fminbnd(@(t) -s.y(t),0,1);
     Beam = s.y(ts);
     Midship = s.x(ts);
     Draft = p.y(fminbnd(p.y,0,1));
     printf('Length: %5.4f \n',LOA)
     printf('Beam: %5.4f \n',Beam)
     printf('Draft: %5.4f \n',Draft)
     printf('Midship: %5.4f \n',Midship)
          
     % Collect everything into a structure
     BezStruct.BezCurves  = BezCurves;
     BezStruct.sections   = sections;
     BezStruct.functions  = functions;
     BezStruct.LOA        = LOA;
     BezStruct.Beam       = Beam;
     BezStruct.Midship    = Midship;
     BezStruct.Draft      = Draft;
    
endfunction

function f = symBezFunc(BezCurves,CrvName)
    % Returns Bezier functions generated using the symbolic package 
    %
    % Inputs:
    %   BezCurves: Structure of Bezier curve control points
    %   CrvName:   Name of curve to generate symbolically 
    % 
    % Output: 
    %    f: Symbolic representation of selected curve
    
    % Points array
    i = cellfind({BezCurves.Curve},CrvName);
    P = [BezCurves(i).x; BezCurves(i).y]';
    
    % Initialize symbolic variable t
    syms t

    % Degree of the curve
    n = size(P,1);

    % Binomial coefficients
    B = binomial(n-1);

    % Terms in t
    t1 = repmat(1-t,n,1);
    t2 = repmat(t,n,1);

    % Powers for each term
    p = (0:n-1)';

    % Raise terms to powers
    t1 = t1 .^ flipud(p);
    t2 = t2 .^ p;

    % Convert Bezier control points to variable precision for symbolic calculations
    P = vpa(P);

    % Symbolic polynomials of Bezier function
    pf = B .* t1 .* t2;
    px = sum(P(:,1) .* pf);
    py = sum(P(:,2) .* pf);

    % Convert to anonymous functions
    f.x = function_handle(px);
    f.y = function_handle(py);
    
endfunction

function t = funcSolve(f,y)    
    % Returns t such that f(t) = y for input values of y                                    
    %                                     
    %                                     
    % Input(s)                            
    %   f: Function handle
    %   y: Values where f(t) are required                                  
    %                                     
    % Output(s)                           
    %   t: Values such that f(t) = y 
    
    % Number of data points
    n = numel(y);
    
    % Initialize output vector
    t = nan(1,n);
    
    % Loop over each y_k to find t_k
    for k = 1:n
        try
            s = fminbnd(@(x) (f(x) - y(k)).^2,0,1);
            t(k) = s;            
        catch
        end_try_catch
    endfor
    
    % Remove values of t outside [0,1]
    t(t < 0) = nan;
    t(t > 1) = nan;

endfunction    


function write2DELFTship(sections,fName)    
    % Input(s)                                  
    %    Hull: Hull points data structure
    %   fName: File name of Excel file, will write output to .txt file
    % Output(s)                                 
    %   Text file with hull cross sections
    
    % Units (feet / meters)
    Use_feet = 1;
    
    % Open file for writing
    delftName = strrep(fName,'xlsx','txt');
    delftName = strrep(delftName,' ','_');
    fid = fopen(delftName,'w');
    
    % Indicate units used
    fprintf(fid,'%d\r\n',Use_feet);
    
    % Space before hull definition points
    fprintf(fid,'\r\n'); 
    
    % Write section curves
    nSections = numel(sections);
    for k = 1:nSections
        x = sections(k).Pts(:,1);
        y = sections(k).Pts(:,2);
        z = sections(k).Pts(:,3);
        for j = 1:numel(x)
            fprintf(fid,'%6.5f  %6.5f  %6.5f\r\n',x(j),y(j),z(j));
        endfor
        % Write blank line
        fprintf(fid,'  \r\n');
    endfor
           
    % Write blank line
    fprintf(fid,'\r\n');    
        
    % Write EOF at end of file
    fprintf(fid,'%s\r\n','EOF');

    fclose(fid);
    
    printf('Section points saved to %s.\n\n', delftName)
    
endfunction    

function z = scale(v,sc)
    % Scales v to values betwen sc(1) and sc(2)
    v0 = min(v(:));
    v1 = max(v(:));
    z  = abs(diff(sc))/(v1 - v0) * (v - v0) + min(sc);
endfunction