matlab对点云进行三角化(las格式读入)

1、LASread_and_convert.m文件（运行这个文件实现las读入为mat文件)

clear;
clc;
close all;
clear all;
A = LASreadAll('house_wall.las');%调用函数读取las文件
Point_City = [A.x,A.y,A.z];
save Point_City.mat;

 %将X1_buny赋值给一个p
    
    p = Point_City;
    
    %显示点云
    figure(1);
hold on
axis equal
title('点云','fontsize',14)
plot3(p(:,1),p(:,2),p(:,3),'g.')%这里颜色如果是 ‘r’就不行？？？
view(-37.5,30)


[t]=MyCrust(p);

%% plot of the oyput triangulation
figure(2)
        hold on
        title('三角化输出','fontsize',14)
        axis equal
        trimesh(t,p(:,1),p(:,2),p(:,3),'edgecolor','b')%plot della superficie trattata
        view(-37.5,30)

2、LASreadAll函数

function A = LASreadAll(infilename)
% LASREADALL reads in all variables from a LAS 1.x data file (used with lidar data) 
%
% INPUT
% infilename: input file name  (for example, 'myinfile.las') 
%           
% OUTPUT
% A:    This is a structure containing all the data in the file.
%         See the file documentation for more information:
% http://www.asprs.org/Committee-General/LASer-LAS-File-Format-Exchange-Activities.html
% 
% EXAMPLE
% A = LASreadAll('infile.las')
%
% Cici Alexander
%     September 2008 (updated 26.09.2008)
% Amy Farris (afarris@usgs.gov) 
%      November 2013 Substatially altered to read in all variables from the file

% This file has the capability to load data files using any 1 of 5 formats.
% However, I was only able to test it with one of the formats. 
% I included support for the rest of the file formats b/c I
% thought someone might find it useful.  I sure hope they work!!!
% If the code crashes, it may be because either the pointDataFormatID 
% value was not what I expected, or the value for pointDataRecordLength 
% is not what I think it should be based on the format ID.  
% If pointDataRecordLength is different than I expected, then you may need 
% to change the values in fseek, eg: (c+????)
% The number added to c should be the sum of bytes of all the variables 
% that occur before the variable currently being read in.
%
% Also, the file I tested this code with was LAS format 1.2, I 
% think this code will run on later versions of LAS.
%
% A brief explanation of the LAS file format:
% LAS files are binary, they begin with header information.  Included in
% the header is the size of the header (in bytes); this is called
% 'OffsetToPointData', refered to a 'c' in this code.  After c bytes the
% data begins with the first x value, then the first y value and so on...
% (exactly what data is included depends on the file format).  Then the
% file continues with the second x value, the second y value and so on...  
% The header tells you how many bytes of data there are for each data point
% ('pointDataRecordLength', also refered to as 'p' in this code).
% So to read in all the x values, you start at the beginig of the file and
% skip c bytes:  "(fseek(fid, c, 'bof');"  
% Then you read in one value and skip p-4 bytes 
% (each x value consists of 4 bytes) and then read the next x and so on:
% "X1 = fread(fid,inf,'int32',p-4);"
%
% This is my (Amy Farris') first attempt at reading in a binary file.  
% I depended strongly on Cici's orignal file (LASRead.m) at first.  Most 
% of the code after about line # 134 was written by me, as are these
% garrolous beginning comments.  I hope they help.

%% Open the file
fid =fopen(infilename);

% Check whether the file is valid
if fid == -1
    error('Error opening file')
end

%% Read in important information from the header
% Check whether the LAS format is 1.1
fseek(fid, 24, 'bof');
VersionMajor = fread(fid,1,'uchar');
VersionMinor = fread(fid,1,'uchar');
% afarris2011Aug20 changed the following line to read LAS1.2 files
% if VersionMajor ~= 1 || VersionMinor ~= 1 
if VersionMajor ~= 1  
    error('LAS format is not 1.*')
end

% Read in the offset to point data
fseek(fid, 96, 'bof');
OffsetToPointData = fread(fid,1,'uint32');

% Read in the point data fotmat ID
fseek(fid, 104, 'bof');
pointDataFormatID = fread(fid,1,'uchar');

% Read in the point data record length
fseek(fid, 105, 'bof');
pointDataRecordLength = fread(fid,1,'short');

% Read in the scale factors and offsets required to calculate the coordinates
fseek(fid, 131, 'bof');
XScaleFactor = fread(fid,1,'double');
YScaleFactor = fread(fid,1,'double');
ZScaleFactor = fread(fid,1,'double');
XOffset = fread(fid,1,'double');
YOffset = fread(fid,1,'double');
ZOffset = fread(fid,1,'double');

% The number of bytes from the beginning of the file to the first point record
% data field is used to access the attributes of the point data
c = OffsetToPointData;

% The number of bytes to skip after reading in each value is based on
% 'pointDataRecordLength' So I need a short version of the variable name:
p = pointDataRecordLength;

%% Now read in the data

% Reads in the X coordinates of the points;  making use of the
% XScaleFactor and XOffset values in the header.
fseek(fid, c, 'bof');
X1 = fread(fid,inf,'int32',p-4);
A.x = X1*XScaleFactor+XOffset;

% Read in the Y coordinates of the points
fseek(fid, c+4, 'bof');
Y1 = fread(fid,inf,'int32',p-4);
A.y = Y1*YScaleFactor+YOffset;

% Read in the Z coordinates of the points
fseek(fid, c+8, 'bof');
Z1 = fread(fid,inf,'int32',p-4);
A.z = Z1*ZScaleFactor+ZOffset;

% Read in the Intensity values of the points
fseek(fid, c+12, 'bof');
A.intensity = fread(fid,inf,'uint16',p-2);

% Read in the Return Number of the points. The first return will have a
% return number of one, the second, two, etc.
% In the next two fread's, p needs to be in bits = 8*byte
fseek(fid, c+14, 'bof');
A.ReturnNumber = fread(fid,inf,'bit3',p*8-3);

% Read in the Number of Returns for a given pulse.
fseek(fid, c+14, 'bof');
fread(fid,1,'bit3');
A.NumberOfReturns = fread(fid,inf,'bit3',p*8-3);

% Read in classification
fseek(fid, c+15, 'bof');
A.classification = fread(fid,inf,'char',p-1);

% Read in scan angle Rank
fseek(fid, c+16, 'bof');
A.scanAngleRank = fread(fid,inf,'char',p-1);
A.scanAngleRankInfo = '-90 to 90 left side';

% Read in User ID
fseek(fid, c+17, 'bof');
A.userID = fread(fid,inf,'char',p-1);

% Read in Point SourceID
fseek(fid, c+18, 'bof');
A.pointSourceID = fread(fid,inf,'short',p-2);


%% Now read in data specific to certain file formats
% Remembeber, only the code that reads in format 1 was tested
switch pointDataFormatID
    case 1 
        if pointDataRecordLength ~= 28
            error('pointDataRecordLength is not what i expected')
        end
        % Read in 'Global Encoding', which tells us what the time variable is
        fseek(fid, 6, 'bof');
        globalEncoding = fread(fid,1,'short');
        if globalEncoding
            A.timeInfo = 'Time is "standard GPS time minus 1e9"';
            disp('Time is "standard GPS time minus 1e9"')
        else
            A.timeIinfo = 'Time is second of the GPS Week; and NO, we do not know WHICH week. ';
            disp('Time is second of the GPS Week; and NO, we do not know WHICH week. ')
        end
        
        % Read in time
        fseek(fid, c+20, 'bof');
        A.time = fread(fid,inf,'double',p-8);
    case 2
        if pointDataRecordLength ~= 26
            error('pointDataRecordLength is not what i expected')
        end
        % Read in color
        fseek(fid, c+20, 'bof');
        A.red = fread(fid,inf,'short',p-2);
        fseek(fid, c+22, 'bof');
        A.green = fread(fid,inf,'short',p-2);
        fseek(fid, c+24, 'bof');
        A.blue = fread(fid,inf,'short',p-2);
    case 3
        if pointDataRecordLength ~= 34
            error('pointDataRecordLength is not what i expected')
        end
        % Read in 'Global Encoding', which tells us what the time variable is
        fseek(fid, 6, 'bof');
        globalEncoding = fread(fid,1,'short');
        if globalEncoding
            A.timeInfo = 'Time is "standard GPS time minus 1e9"';
            disp('Time is "standard GPS time minus 1e9"')
        else
            A.timeInfo = 'Time is second of the GPS Week; and NO, we do not know WHICH week. ';
            disp('Time is second of the GPS Week; and NO, we do not know WHICH week. ')
        end
        
        % Read in time
        fseek(fid, c+20, 'bof');
        A.time = fread(fid,inf,'double',p-8);
        
        % Read in color
        fseek(fid, c+28, 'bof');
        A.red = fread(fid,inf,'short',p-2);
        fseek(fid, c+30, 'bof');
        A.green = fread(fid,inf,'short',p-2);
        fseek(fid, c+32, 'bof');
        A.blue = fread(fid,inf,'short',p-2);
    case 4
        if pointDataRecordLength ~= 57
            error('pointDataRecordLength is not what i expected')
        end
        % Read in 'Global Encoding', which tells us what the time variable is
        fseek(fid, 6, 'bof');
        globalEncoding = fread(fid,1,'short');
        if globalEncoding
            A.timeInfo = 'Time is "standard GPS time minus 1e9"';
            disp('Time is "standard GPS time minus 1e9"')
        else
            A.timeInfo = 'Time is second of the GPS Week; and NO, we do not know WHICH week. ';
            disp('Time is second of the GPS Week; and NO, we do not know WHICH week. ')
        end
        
        % Read in time
        fseek(fid, c+20, 'bof');
        A.time = fread(fid,inf,'double',p-8);
        
        fseek(fid, c+28, 'bof');
        A.wavePacketDescriptorIndex = fread(fid,inf,'char',p-1);
        fseek(fid, c+29, 'bof');
        % i am not sure that 'unit64' is correct below
        A.byteOffsettoWaveformData = fread(fid,inf,'unit64',p-8);
        fseek(fid, c+37, 'bof');
        A.waveformPacketSize = fread(fid,inf,'unit32',p-4);
        A.waveformPacketSizeInfo = 'in bytes';
        fseek(fid, c+41, 'bof');
        % i am not sure that 'float' is correct below
        A.returnPointWaveformLocation = fread(fid,inf,'float',p-4);
        fseek(fid, c+45, 'bof');
        A.Xt = fread(fid,inf,'float',p-4);
        fseek(fid, c+49, 'bof');
        A.Yt = fread(fid,inf,'float',p-4);
        fseek(fid, c+53, 'bof');
        A.Zt = fread(fid,inf,'float',p-4);
    case 5
        if pointDataRecordLength ~= 63
            error('pointDataRecordLength is not what i expected')
        end
        % Read in 'Global Encoding', which tells us what the time variable is
        fseek(fid, 6, 'bof');
        globalEncoding = fread(fid,1,'short');
        if globalEncoding
            A.timeInfo = 'Time is "standard GPS time minus 1e9"';
            disp('Time is "standard GPS time minus 1e9"')
        else
            A.timeInfo = 'Time is second of the GPS Week; and NO, we do not know WHICH week. ';
            disp('Time is second of the GPS Week; and NO, we do not know WHICH week. ')
        end
        
        % Read in time
        fseek(fid, c+20, 'bof');
        A.time = fread(fid,inf,'double',p-8);
        
        % Read in color
        fseek(fid, c+28, 'bof');
        A.red = fread(fid,inf,'short',p-2);
        fseek(fid, c+30, 'bof');
        A.green = fread(fid,inf,'short',p-2);
        fseek(fid, c+32, 'bof');
        A.blue = fread(fid,inf,'short',p-2);
        fseek(fid, c+34, 'bof');
        A.wavePacketDescriptorIndex = fread(fid,inf,'char',p-1);
        fseek(fid, c+35, 'bof');
        % i am not sure that 'unit64' is correct below
        A.byteOffsettoWaveformData = fread(fid,inf,'unit64',55);
        fseek(fid, c+43, 'bof');
        A.waveformPacketSize = fread(fid,inf,'unit32',p-4);
        A.waveformPacketSizeInfo = 'in bytes';
        fseek(fid, c+47, 'bof');
        % i am not sure that 'float' is correct below
        A.returnPointWaveformLocation = fread(fid,inf,'float',p-4);
        fseek(fid, c+51, 'bof');
        A.Xt = fread(fid,inf,'float',p-4);
        fseek(fid, c+53, 'bof');
        A.Yt = fread(fid,inf,'float',p-4);
        fseek(fid, c+57, 'bof');
        A.Zt = fread(fid,inf,'float',p-4);
end

3、MyCrust函数

%% MyCrust
%
%Simple surface recostruction program based on Crust algorithm
%Given a set of 3D points returns a triangulated tight surface.
%
%The more points there are the best the surface will be fitted,
%although you will have to wait more. For very large models an 
%help memory errors may occurs.
%It is important even the point distribution, generally uniformly 
% distributed points with denser zones in high curvature features
% give the best results.
% 
% Remember  crust algorithom needs a cloud representing a volume 
% so open surface may give inaccurate results. 
%  
%  
% If any problems occurs in execution, or if you found a bug,
% have a suggestion or question just contact me at:  
%   
% giaccariluigi@msn.com
% 
%
%
%
%Here is a simple example:
%
%load Dino.mat%load input points from mat file
%
%[t]=MyCrust(p);
%
% figure(1)
%         hold on
%         title('Output Triangulation','fontsize',14)
%         axis equal
%         trisurf(t,p(:,1),p(:,2),p(:,3),'facecolor','c','edgecolor','b')
%
%Input:
%              p is a Nx3 array containing the 3D set of points
%Output:
%              t are points id contained in triangles Nx3 array too.
%
% See also qhull, voronoin, convhulln, delaunay, delaunay3, tetramesh.
%
%Author:Giaccari Luigi
%Last Update: 1/12/2008
%Creation: 10/10/2008


function [t]=MyCrust(p)
%%   Main
starttime=clock;


%add points to the given ones, this is usefull
%to create outside tetraedrom
tic
p=AddShield(p);
fprintf('Addedded Shield: %4.4f s\n',toc)


tic
tetr=delaunayn(p);%creating tedraedron
tetr=int32(tetr);%use integer to save memory
fprintf('Delaunay Triangulation Time: %4.4f s\n',toc)


%connectivity data
%find triangles to tetraedrom and tetraedrom to triangles connectivity data
tic
[t2tetr,tetr2t]=Connectivity(tetr);
fprintf('Connectivity Time: %4.4f s\n',toc)
tic
[cc,r]=CC();%Circumcenters of tetraedroms
fprintf('Circumcenters Time: %4.4f s\n',toc)
clear n



tic
t=Walking();%Flagging tetraedroms as inside or outside
fprintf('Walking Time: %4.4f s\n',toc)


time=etime(clock,starttime);
fprintf('Total Time: %4.4f s\n',time)





%% Circumcenters(Nested)
    function [cc,r]=CC()
         %finds circumcenters fro a set of tetraedrom
        
         %points of tetraedrom 
        p1=(p(tetr(:,1),:));
        p2=(p(tetr(:,2),:));
        p3=(p(tetr(:,3),:));
        p4=(p(tetr(:,4),:));

        %vectors of tetraedrom edges
        v21=p(tetr(:,1),:)-p(tetr(:,2),:);
        v31=p(tetr(:,3),:)-p(tetr(:,1),:);
        v41=p(tetr(:,4),:)-p(tetr(:,1),:);

        %preallocation
        cc=zeros(size(tetr,1),3);






         %Solve the system using cramer method
        d1=sum(v41.*(p1+p4)*.5,2);
        d2=sum(v21.*(p1+p2)*.5,2);
        d3=sum(v31.*(p1+p3)*.5,2);

        det23=(v21(:,2).*v31(:,3))-(v21(:,3).*v31(:,2));
        det13=(v21(:,3).*v31(:,1))-(v21(:,1).*v31(:,3));
        det12=(v21(:,1).*v31(:,2))-(v21(:,2).*v31(:,1));

        Det=v41(:,1).*det23+v41(:,2).*det13+v41(:,3).*det12;



        detx=d1.*det23+...
            v41(:,2).*(-(d2.*v31(:,3))+(v21(:,3).*d3))+...
            v41(:,3).*((d2.*v31(:,2))-(v21(:,2).*d3));



        dety=v41(:,1).*((d2.*v31(:,3))-(v21(:,3).*d3))+...
            d1.*det13+...
            v41(:,3).*((d3.*v21(:,1))-(v31(:,1).*d2));



        detz=v41(:,1).*((v21(:,2).*d3)-(d2.*v31(:,2)))...
            +v41(:,2).*(d2.*v31(:,1)-v21(:,1).*d3)...
            +d1.*(det12);

        %Circumcenters
        cc(:,1)=detx./Det;
        cc(:,2)=dety./Det;
        cc(:,3)=detz./Det;


         %Circumradius
        r=((sum((p2-cc).^2,2))).^.5;%quadrato del raggio
    end



%% Connectivity(Nested)

    function [t2tetr,tetr2t]=Connectivity(tetr)

       

        numt = size(tetr,1);
        vect = 1:numt;
        t = [tetr(:,[1,2,3]); tetr(:,[2,3,4]); tetr(:,[1,3,4]);tetr(:,[1,2,4])];%triangles not unique
        [t,j,j] = unique(sort(t,2),'rows');%triangles
        t2tetr = [j(vect), j(vect+numt), j(vect+2*numt),j(vect+3*numt)];%each tetraedrom has 3 triangles


        % triang-to-tetr connectivity
      
        nume = size(t,1);
        tetr2t  = zeros(nume,2,'int32');
        count= ones(nume,1,'int8');
        for k = 1:numt

            for j=1:4
                ce = t2tetr(k,j);
                tetr2t(ce,count(ce)) = k;
                count(ce)=count(ce)+1;
            end

        end


    end      % connectivity()



%% Walking(Nested)
    function t=Walking()
        np=size(p,1)-540;%540 = number of shield points put at the end of array
        numtetr=size(tetr,1);
        nt=size(tetr2t,1);


        deleted=true(numtetr,1);%deleted tetraedroms
        checked=false(numtetr,1);%checked tetraedros 
        onfront=false(nt,1);%tetraedroms that need to be checked
        countchecked=0;
        %First flag as outsde tetraedroms with Shield points
        for i=1:numtetr
            for j=1:4
                if tetr(i,j)>np;
                    deleted(i)=true;
                    checked(i)=true;
                    onfront(t2tetr(i,:))=true;
                    countchecked=countchecked+1;
                    break
                end
            end
        end

        %tollerances to mark as in or out
        toll=zeros(nt,1)+.95; 
        level=0;
        
        alpha=zeros(nt,1);%intersection factor
        %it is computed from radius of the tetraedroms circumscribed sphere
        % and the distance between their center
        for i=1:nt
            if tetr2t(i,2)>0 %jump boundary tetraedrom
                    distcc=sum((cc(tetr2t(i,1),:)-cc(tetr2t(i,2),:)).^2,2);%distance from circumcenters
                    %intersection factor
                    alpha(i)=(-distcc+r(tetr2t(i,1))^2+r(tetr2t(i,2))^2)/(2*r(tetr2t(i,1))*r(tetr2t(i,2)));
            end
        end
        clear cc
        
%         Now we scan all tetraedroms. When one is scanned put on front is
%         neighobur. This means that now even the neighobour can be
%         checked. At the begining only tetraedroms with shield points are
%         on front, because we are sure the are out.
%         Tetraedrom with high intersction factor  will be marked as equal
%         else different.
%         When I say high i mean under a set tollerance that becames lower as the algorithm
%         progresses. This Aims to avoid errors propagation when a tetraedrom is wrong marked.
%         
        
        while countchecked<numtetr && level<50
             level=level+1;%level of scan reached
             
            for id=1:nt%loop trough triangles

                if onfront(id)

                    tetr1=tetr2t(id,1);tetr2=tetr2t(id,2);%tetraedroms linked to triangle under analysis
                    if  tetr2==0 %do not check boundary triangles
                        onfront(id)=false;
                        continue
                       
                    elseif (checked(tetr1) && checked(tetr2)) %tetraedroms are already checked
                         onfront(id)=false;
                        continue
                       
                    end

                    if alpha(id)>=toll(id) %flag as equal
                        if checked(tetr1)%find the checked one between the two
                            deleted(tetr2)=deleted(tetr1) ;%flag as equal
                            checked(tetr2)=true;%check
                            countchecked=countchecked+1;
                            onfront(t2tetr(tetr2,:))=true;%put on front all tetreadrom triangles
                        else
                            deleted(tetr1)=deleted(tetr2) ;%flag as equal
                            checked(tetr1)=true;%check
                            countchecked=countchecked+1;
                            onfront(t2tetr(tetr1,:))=true;%put on front all tetreadrom triangles
                        end
                         onfront(id)=false;%remove from front
                   
                    
                    elseif alpha(id)<-toll(id)%flag as different
                        if checked(tetr1)%find the checked one between the two
                            deleted(tetr2)=~(deleted(tetr1)) ;%flag as different
                            checked(tetr2)=true;%check
                            countchecked=countchecked+1;
                            onfront(t2tetr(tetr2,:))=true;%put on front all tetreadrom triangles
                        else
                            deleted(tetr1)=~(deleted(tetr2)) ;%flag as different
                            checked(tetr1)=true;%check
                            countchecked=countchecked+1;
                            onfront(t2tetr(tetr1,:))=true;%put on front all tetreadrom triangles
                        end
                         onfront(id)=false;%remove from front

                         
                    else
                        toll(id)=toll(id)-.05;%tolleraces were too high next time will be lower

                    end



                end
            end
            
            if level==31 %brute continuation(this may appens when the triangulation is corrupt)
                 warning('Brute continuation necessary')
                onfront(t2tetr(~(checked),:))=true;%force onfront collocation
            end
        end
        
        %delete tetraedroms marked as outside 
        tetr(deleted,:)=[];
        
        
        
        %take boundary triangles from tetraedroms
        %this is the raw surface and needs improvements to be used in CAD
        %systems. Maybe in my next revision I will add surface post treatments.
        %Anyway for grafical purpose this should be good.

        %extract boundary triangles 
        t=BoundTriangles(tetr);
        
        
        %Output Data
        numchecked=countchecked/numtetr;
        if level==50
            warning([num2str(level),' th level was reached\n'])
        else
         fprintf('%4.4f th level was reached\n',level)
        end
        fprintf('%4.4f %% of Tetraedrom were checked\n',numchecked*100)



    end




end

%% AddAddShield
function pnew=AddShield(p)

%find the bounding box
maxx=max(p(:,1));
maxy=max(p(:,2));
maxz=max(p(:,3));
minx=min(p(:,1));
miny=min(p(:,2));
minz=min(p(:,3));

%give offset to the bounding box
step=max(abs([maxx,minx,maxy,miny,maxz,minz]));

maxx=maxx+step;
maxy=maxy+step;
maxz=maxz+step;
minx=minx-step;
miny=miny-step;
minz=minz-step;
step=step/100;
N=10;



%creating a grid lying on the bounding box
vx=linspace(minx,maxx,N);
vy=linspace(miny,maxy,N);
vz=linspace(minz,maxz,N);




[x,y]=meshgrid(vx,vy);
facez1=[x(:),y(:),ones(N*N,1)*maxz];
facez2=[x(:),y(:),ones(N*N,1)*minz];
[x,y]=meshgrid(vy,vz-step);
facex1=[ones(N*N,1)*maxx,x(:),y(:)];
facex2=[ones(N*N,1)*minx,x(:),y(:)];
[x,y]=meshgrid(vx-step,vz);
facey1=[x(:),ones(N*N,1)*maxy,y(:)];
facey2=[x(:),ones(N*N,1)*miny,y(:)];

%add points to the p array
pnew=[p;
 facex1;
 facex2;
 facey1;
 facey2;
 facez1;
 facez2];

% figure(2)
% plot3(pnew(:,1),pnew(:,2),pnew(:,3),'.g')


end



%% BoundTriangles
function t=BoundTriangles(tetr)

numt = size(tetr,1);
vect = 1:numt;
t = [tetr(:,[1,2,3]); tetr(:,[2,3,4]); tetr(:,[1,3,4]);tetr(:,[1,2,4])];
[t,j,j] = unique(sort(t,2),'rows');
t2tetr = [j(vect), j(vect+numt), j(vect+2*numt),j(vect+3*numt)];


% triang-to-tetr connectivity
% Each row has two entries corresponding to the triangle numbers
% associated with each triangle. Boundary triangles have only one tetraedrom
nume = size(t,1);
tetr2t  = zeros(nume,2,'int32');
count= ones(nume,1,'int32');
for k = 1:numt

    for j=1:4
        ce = t2tetr(k,j);
        tetr2t(ce,count(ce)) = k;
        count(ce)=count(ce)+1;
    end

end
tbound=false(numt,1);
tbound=tetr2t(:,2)==0;
t=t(tbound,:);
end