function retval = FEMKonstructor(code,subcode,data)
% 
% Let's deal with two objects: One for the interface.
% 
% A new page is opened: The old version of this chaotic monster is now oldFEMKonstructor.m
% I will now screw through this hack and first of all delete all functions that
% are hideous and need replacement. 
%
if nargin<3, data=[]; end;
if nargin<2, subcode=[]; end;

windowobj = findobj('tag','FEMKonstructor');
if isempty(windowobj), Kmessage('FEMKonstructor could not find data structure'); return; end;
Konst = get(windowobj,'userdata');
savebackinterface = 1;

geometryhandle = findobj('tag','FEMGeometryobject');
if ~isempty(geometryhandle),
  mobj = get(geometryhandle,'userdata');
else
  mobj=makeGeometry;
  Kmessage('make new geometry object');
  geometryhandle = makestoreobject('FEMGeometryobject',mobj,'disp(get(gcbo,''userdata''));');
end;

mobj.dosave=1;

if Konst.axes ~= gca, 
  delete(gca);                
  axes(Konst.axes); 
end;    % switch always.

%============================ 'switch' ===========================

switch code,
 case 'cludge',
  mobj=make_a_cludge(Konst,mobj);   % Why would this thing be first ?? .;-)

 case 'addvertexlist',
  mobj=add_vertex_list(Konst,mobj,subcode);
  
 case 'record',
  Kmessage_recorder(subcode);
  
 case 'mouse',    %this is called si rien d'autre a passé.
  mobj=mouse_operation(Konst,mobj);
  
 case 'setmousecode', 
  mobj.mouseoperation = subcode;
  
 case 'gluonset',
  mobj.gluons = subcode;    % 0 or 1: postpone glueing by storing info in a dotted line.
  
 case 'gluonshow',
  show_gluons(mobj,subcode,data);
  
 case 'top',    axes(Konst.axes); view([0 90]);
 case 'bottom', axes(Konst.axes); view([0 -90]);
 case 'front',  axes(Konst.axes); view([0 0]);
 case 'back',   axes(Konst.axes); view([180 0 ]);
 case 'left',   axes(Konst.axes); view([90 0]);  
 case 'right',  axes(Konst.axes); view([-90 0]);  
 case 'rotate', axes(Konst.axes); rotate3d;
 case 'projection', 
  switch subcode,
   case 0, set(Konst.axes,'projection','orthographic');
   case 1, set(Konst.axes,'projection','perspective');
  end;
  
 case 'zoom',   axes(Konst.axes); mobj=zoominout(Konst,mobj,subcode);
 case 'viewsel',
  Konst = select_view(Konst,mobj,subcode);
  
 case 'setaxisrange',  axes(Konst.axes); axis('equal'); 
  Konst.axisrange=axis; 
  xlabel('X'); ylabel('Y'); zlabel('Z');
  set(windowobj,'userdata',Konst);
  
 case 'color',
  mobj=change_curve_color(Konst,mobj,subcode);
  
 case 'xyz',    % something happened with the xyz panel.
  mobj= change_xyz_values(Konst,mobj);
  if mobj.autoupdate>0,
    mobj = update_hierarchy(Konst,mobj);
  end;
 case 'projectpoint',
  mobj= project_point_inplane(Konst,mobj);

 case 'projectcurve',
  mobj= project_curve_inplane(Konst,mobj);
  
 case 'autoproject',
  mobj = toggle_autoproject(Konst,mobj,subcode);   % automatic projection on active plane.

 case 'nearplane',
  mobj = activate_near_points(Konst,mobj,subcode);
  
 case 'update',
  mobj = update_hierarchy(Konst,mobj);    % reaction to key 'u' or 'U'
  
 case 'geometry',
  Kmessage('FEMKonstructor geometry');
  
 case 'load',
  [mobj,flg] = load_a_file(Konst,mobj,subcode);
  
 case 'points',        mobj=point_indexing(Konst,mobj,subcode); 
 case 'curves',        mobj=curve_indexing(Konst,mobj,subcode);
 case 'surfaces',      Kmessage([code,' ', subcode]);
 case 'cells',         Kmessage([code,' ', subcode]);

 case 'showsurf',
  mobj=make_surfaces(Konst,mobj);

 case 'showskin',
  mobj=make_skinsurfaces(Konst,mobj);
  
 case 'killsurf',
  kill_surfaces(mobj);
  
 case 'light',
  if subcode==0,cod='off'; else,cod='on';end;
  mobj=light_operation(mobj,cod);
  
 case 'makelight',  
  mobj=light_operation(mobj,'on');
  
 case 'killlight',
  mobj=light_operation(mobj,'off');
 
 case 'lightpos',
  mobj=light_operation(mobj,'move',subcode);
  
 case 'beginspline',
  mobj=begin_spline(Konst,mobj);
  
 case 'killcurve',
  mobj=destroy_current_curve(Konst,mobj);

 case 'killpoint',
  mobj=kill_point(Konst,mobj);

 case 'memorize',
  mobj = memorize_obj(Konst,mobj,subcode,data);
  
 case 'recallpoints',
  retval=recall_points(Konst,mobj);

 case 'showmemory',
  retval = get_memory_content(mobj);

 case 'forget',
  mobj = clear_object_memory(Konst,mobj,subcode); % subcode 'points','curves','loops'
  
 case 'showmemory',
  mobj = show_memory(Konst,mobj,subcode);  % subcode 0 or 1
  
 case 'showtopology',
  mobj = show_vertices_numbers(Konst,mobj,subcode);  % subcode 0 or 1
 
 case 'cellnumbers',    % call 'by hand'
  mobj = show_cell_numbers(Konst,mobj,subcode,data);
  
 case 'cellno',         % call 'by button callback'
  if subcode==1,
    mobj = show_cell_numbers(Konst,mobj);
  elseif subcode==0,
    mobj = show_cell_numbers(Konst,mobj,'off');
  end;
  
% case 'cellvisible',
%  mobj = cells_visibles(Konst,mobj,subdata,data);
  
 case 'showsplinedir',
  mobj = show_spline_direction(Konst,mobj,subcode);
  
 case 'attach',
  [mobj,flg] = attach_pt_to_spline(Konst,mobj,subcode); % subcode 1: at end, subcode 0: at beginning
  
 case 'joinsplines',
  mobj = join_connected_splines(Konst,mobj);
  
 case 'flipmemory',
  mobj = permute_memory(Konst,mobj);    % uses edit field tagged editmemorysequence'
  
 case 'makevertexpoint', % turn current redpoint into vertex.
  mobj = turn_into_vertex(Konst,mobj);
  
 case 'makecurvepoint',
  mobj = make_curve_point(Konst,mobj);
  
 case 'selectable', 
  mobj=set_selectables(Konst,mobj,subcode,data);  

 case 'preselect',
  mobj=set_preselect(Konst,mobj,subcode);  % subcode: either 'points' or 'curves';
  
 case 'hit',
  mobj = hit_object(Konst,mobj,subcode,data);  % clicked on point or curve.
  
 case 'make',
  mobj = select_point(Konst,mobj,0);
  xyz=getedit_values(Konst.xyzedits);
  [mobj,flg]= addobject(Konst,mobj,'point',xyz,1);
  mobj = select_point(Konst,mobj,flg);
  
 case 'makecpoint',
  [mobj,flg]= insert_curve_point(Konst,mobj);
  
 case 'splitspline',  
  [mobj,flg]= split_curve_at_pt(Konst,mobj);   % 
  
 case 'midpoint',
  mobj = newpoint_between_two(Konst,mobj,subcode);  % subcode contains number of intermediate points.

 case 'gluepoint',
  if isempty(subcode),
    [mobj,flg]= glue_point_to_point(Konst,mobj);
  else
    [mobj,flg]= glue_point_to_point(Konst,mobj,subcode(1),subcode(2));
  end;    
  
 case 'glueblock',
  if mobj.gluons==1,
    mobj = glue_gluon_pairs(Konst,mobj);
  else
    [mobj,flg]=glue_near_blocks(Konst, mobj); 
  end;
  
 case 'rerender',
  mobj= rerender_graph(Konst,mobj);

 case 'reduce',  
  if subcode==0,cod='off',else,cod='on';end;
  minimize_visible(mobj,cod);
  
 case 'createtopology',
  name = get(Konst.namefield,'string');
  name = make_a_label(name,8);
  [mobj,flg]= create_free_cell(Konst,mobj,subcode,name);
  
 case 'bodymodify',
  [mobj,flg]= move_free_cell(Konst,mobj);
  
 case 'create',
  mobj = create_object(Konst,mobj,subcode);

 case 'startobject',
  mobj = start_object(Konst,mobj);    % such as spline, cube etc.  

 case 'remake',
  mobj = load_picture(Konst,subcode,0);        % start from scratch

 case 'reload',
  mobj = load_picture(Konst,subcode,1);        % add to existing
  
 case 'save', 
  mobj = save_picture(Konst,mobj,subcode);   % subcode is filename

 case 'savemacro', 
  mobj = save_macro(Konst,mobj,subcode);   % subcode is filename

 case 'loadmacro', 
  [mobj,retval] = include_macro(Konst,mobj,subcode,data);   % subcode is filename

 case 'savecurrentmac', 
  subcode=nameregistration('new');
  mobj = save_macro(Konst,mobj,subcode);   % subcode is filename

 case 'loadcurrentmac', 
  subcode=nameregistration('last');
  [mobj,retval] = include_macro(Konst,mobj,subcode,[]);   % subcode is filename

 case 'savebin',
  mobj = save_binary(Konst,mobj);
  
 case 'loadbin',
  mobj = load_binary(Konst);   

 case 'saveselect', 
  mobj = save_selected_pict(Konst,mobj,subcode);   % subcode is filename
  
 case 'add',
  [mobj,flg]= addobject(Konst,mobj,subcode,data);
  
 case 'glue',
  [mobj,flg]= show_topology(Konst,mobj,subcode);   % turn memorypoints into tetra, prism or cube
  
 case 'popoutcube',
  mobj = popout_cube(Konst,mobj);

 case 'popoutprism',
  if subcode==3,
  mobj = popout_prism(Konst,mobj);
  elseif subcode==4,
  mobj = popout_roof(Konst,mobj);
  end;
  
  
 case 'tissuename',
  set_tissue_namefield(Konst,subcode);
  
 case 'translateall',
  mobj=translate_all(Konst,mobj,subcode);
  
 case 'rotateall',
  mobj=rotate_all(Konst,mobj,subcode);
  
 case 'scaleall',
  mobj=scale_all(Konst,mobj,subcode);

  % Obsolete stuff, but may be re-used later.    
 case 'bezier',
  Kmessage('convert to bezier');
  mobj = convert_to_bezier(Konst,mobj);
  Kmessage('render bezier');
  mobj = render_all_bezier(Konst,mobj);
  Kmessage('make edges');
  mobj = make_bezedges(Konst,mobj);
  toggle_oldlines('off');
  
 case 'bezierselect',
  mobj = bezier_select(Konst,mobj,subcode);
  
 case 'seespline',
  toggle_oldlines(subcode);
  
 otherwise,
  disp('FEMKonstructor: Illegal code');
  
end;  
%============================ switch over  ===========================
if savebackinterface==1,
  set(windowobj,'userdata',Konst);
end;

if mobj.dosave==1,
  set(geometryhandle,'userdata',mobj);
end;
return;  % end of main program.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function  M=mouse_operation(Ko,M)
  but = analyse_mouse_event(Ko.fig);   % but is button number
  switch M.mouseoperation,
   case 0,           % nothing
    ret=ebenenmaster('clicked');
    check=ret{1};
    if check==1,
      xyz=ret{2};
      xyz= round(xyz*10)/10;
      Kmessage(sprintf('%6.1f %6.1f %6.1f',xyz(1:3)),1);
    else,
      Kmessage('');
    end;
   
   case 1,
    if but==2, M = move_to_projection(Ko,M,0); end; % move currently selected point to mouse point (on surface).
   
   case 2, 
    if but==2, M = move_to_projection(Ko,M,1); end; % move selected point parallel to
                                                    % surface
   case 3,
    ret=ebenenmaster('clicked');
    check=ret{1};
    if check==1,
      xyz=ret{2};
      xyz= round(xyz*10)/10;
      mss = datacollection(but,xyz);        % external for now.
      if ~isempty(mss),
	Kmessage(mss);
      end;
    end;
    
   case 4,  % fiber direction
    ret = ebenenmaster('clicked');
    check=ret{1};
    if check==1,
      xyz= ret{2};
      [iii,dire,sif] = fiberdirection(xyz,8,1);
      dire=dire';
      if ~isempty(iii) & iii(1) > 0,
	%     Kmessage(sprintf('coords: %5.2f %5.2f %5.2f .  Indicies %4i %4i %4i . Dir %5.2f
	%     %5.2f  %5.2f s: %6.3f',xyz,iii,dire,sif));
	Kmessage(sprintf('Significance: %6.3f',sif));
	o = findobj('tag','muscledirection');
	if ~isempty(o) & ishandle(o), delete(o); end;
	o=[0 0 0];
	xdir = [xyz(:)'-dire(:,1)';xyz(:)'+dire(:,1)'];
	o(1) = line('xdata',xdir(:,1),'ydata',xdir(:,2),'zdata',xdir(:,3),'color','y');
	xdir = [xyz(:)'-dire(:,2)';xyz(:)'+dire(:,2)'];
	o(2) = line('xdata',xdir(:,1),'ydata',xdir(:,2),'zdata',xdir(:,3),'color','g');
	xdir = [xyz(:)'-dire(:,3)';xyz(:)'+dire(:,3)'];
	o(3) = line('xdata',xdir(:,1),'ydata',xdir(:,2),'zdata',xdir(:,3),'color','m');
	set(o,'tag','muscledirection');
      end;
    end;
   case 5,     % generate a vertex at current position.
    M = select_point(Ko,M,0);
    xyz=get_projection_coords(Ko,M);
    if ~isempty(xyz),
      [M,np ]= addobject(Ko,M,'point',xyz,1);
      M = select_point(Ko,M,np);    
      M.vertextable = [M.vertextable,np];
    end;
   case 6,     % add a point to a data structure that is associated with
               % the current ebene.
       ret=ebenenmaster('clicked');
       check=ret{1};
       if but==2,
	 fprintf(1,'S.C0=[ %7.3f %7.3f %7.3f]; \n',ret{5});
	 fprintf(1,'S.V1=[ %7.3f %7.3f %7.3f]; \n',ret{6});
	 fprintf(1,'S.V2=[ %7.3f %7.3f %7.3f]; \n',ret{7});
       elseif but==3,
	 fprintf(1,'\n S.N=0;\n');
       elseif but==1,
	 if check==1,
	   fprintf(1,'S.N=S.N+1; S.V(S.N,1:2)=[%7.3f %7.3f ];\n',ret{3}(1),ret{3}(2));
	 else 
	   fprintf(1,'\n ');
	 end;
       end;
  end;
  
  return;

%----------------------------------  %
  
function M=add_vertex_list(Ko,M,v)
  [ll,n3]=size(v);
  for k=1:ll,
    [M,np]=addobject(Ko,M,'point',v(k,1:3),1);
    M = select_point(Ko,M,np);
    M.vertextable = [M.vertextable,np];
  end;
  return;
  

%-----------------------------------------------------------
function but=analyse_mouse_event(fig)
x=get(fig,'selectiontype');
switch x,
case 'normal', but=1;
case 'extend', but=2;
case 'alt',    but=3;
case 'open',   but=4;
end;
return;
%-------------------------------------------------

function xyz=getedit_values(edits)
   for k=1:length(edits),
       xyz(k) = str2num(get(edits(k),'string'));
   end;
return;

%-----------------------------------------------------------
function setedit_values(edits,xyz)
   for k=1:length(edits),
       set(edits(k),'string',sprintf('%4.1f',xyz(k)));
   end;
return;

%-----------------------------------------------------------
% selection of curves: If a new curve is selected, then point selection 
%should mean only point selection of the curve points. 
%same with surface points. 
%Interior points?

%Wath does geometrix do that FEMKonstructor does not? 
% Answer: geometrix was totally ignored just like Troubadix.
%-----------------------------------------------------------

function o =makeGeometry
o.dosave=1;
o.autoupdate=1;
o.freecell=0;         % gnadenlos....
o.freecellnumber=0;   % impitoyablement ...
o.specialcell=0;
o.macrocoords=[];     % merciless...  resetreset.
o.mouseoperation=0;
o.gluons=0;
o.bezierdefined=0;
o.pointhandles=[];
o.curvehandles=[];
o.surfhandles=[];
o.vertextable=[];
o.edgetable=[];
o.edgehandles=[];
o.pointtypes=[];
o.pointreferences=[];    % NOTE there is something seriously wrong with this, should be a curve number in it for each point.
o.curvetypes=[];
o.curvereferences=[];
o.xyz = [];
o.A = sparse([]);
o.B = sparse([]);
o.C = sparse([]);
o.selectablepoints=[]; 
o.selectedpoints=[];
o.selectedcurves=[];
o.redpoint=0;
o.lastpoint=0;
o.redcurve=0;
o.lastcurve=0;
o.redsurface=0;
o.redcell=0;
o.selectedptindex=0;
o.curves={};        % qcurve handles
o.loopintern={};
o.cells={};
o.celllabels={};
o.memory=[];
o.memoryperm=[];
o.curvememory=[];
o.loopmemory=[];
o.splinepointindex=0;
return;

%-----------------------------------------------------------
function [M,flg]= addobject(Ko,M,co,data,typus,reference)
% create points or edges (old type) co is either 'point' or 'curve'. 
% attention, this method was partially abandoned in further developments.
% but the last thing that I added was the generation of Bezier points on edges and
% surfaces. For them the types are 3 and 4 and the reference field referes to edge number
% or surface number.
% 
if nargin<5, typus=0; end;
if nargin<6, reference=0; end;

flg=1;

switch co,
 case 'point',     % data is xyz coordinates
     [msize,colo,mark]=objproperties('point',typus,0);
     h = ...
      line('xdata',data(1),'ydata',data(2),'zdata',data(3),'marker',mark,'markersize',...
     msize,'linestyle','none','color',colo);
     M.pointhandles = [M.pointhandles;h];
     M.pointtypes = [M.pointtypes;typus];
     M.pointreferences = [M.pointreferences;reference];   % reference to edge or surface or whatever
     switch typus,
        case 0,  adde = [ data(1),data(2),data(3), zeros(1,3)];
	case 1,  adde = [ data(1),data(2),data(3), zeros(1,3)];
	case 2,  adde = [ data(1),data(2),data(3), data(4), 0 , 0];
	case 3,  adde = [ data(1),data(2),data(3), 0 , 0 , 0];  % control or Bezier control points on an edge
	case 4,  adde = [ data(1),data(2),data(3), 0 , 0 , 0];  % Bezier points on surface
     end;
     M.xyz = [M.xyz; adde];
     np = length(M.pointhandles);     flg=np;          % pointindex.
     set(h,'tag','point','userdata',[np,typus]);
     flg=np;

 case 'curve',     % data contains point numbers. 
      % get xyz      for k=1:length(data),
      for k=1:length(data),
	np = data(k);
	if np>0 & np <= length(M.pointhandles),
	  h = M.pointhandles(np);
	  xyz(k,:) = [get(h,'xdata'),get(h,'ydata'),get(h,'zdata')];
	else
	  Kmessage('add curve: illegal point number');
	  flg=0;
	end;  
      end;
      xyz = cubicsplinecurve(xyz,length(data)*5,'equal');
      [msize,colo,mark]=objproperties('curve',typus,0);
      h = line('xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3),'linewidth',msize,'color',colo);
      M.curvehandles = [M.curvehandles;h];
      M.curvetypes=[M.curvetypes;typus];
      M.curvereferences=[M.curvereferences;reference];
      ncur = length(M.curvehandles);
      flg=ncur;
      set(h,'userdata',ncur,'tag','curve');
%      disp(length(M.pointhandles));
      
      nc=length(M.curves);
      M.curves{nc+1} = data;
      M=update_curvepoint(Ko,M,nc+1);        % if there are curve points on the curve, update them. 

      % change the references field for the points
      for j=1:length(data),
	ix = data(j);
	ty=M.pointtypes(ix);
	if ty~=2,        % curve points
	  M.pointreferences(ix) = 1;  % no longer free vertex.
	end;
      end;
    end;  % case
     
return;
    
%-----------------------------------------------------------
function [M,ret]= load_a_file(Ko,M,filename);

fid = fopen(filename,'r');
ret=1;
eof=0;

nptprior = length(M.pointhandles);
while ~eof,
  [C,n] = fscanf(fid,'%c',1);
  switch C,
  case 'E',
  eof=1;
  break;
  case 'P',
%  [no,n]=fscanf(fid,'%i',1);
  [xyz,n]=fscanf(fid,'%f',3);
  [M,ret]= addobject(Ko,M,'point',xyz,0);
%  disp(xyz');

  case 'C',
  [no,n] = fscanf(fid,'%i',1);
  [no,n] = fscanf(fid,'%i',no);
%  disp(no');
  no = no+nptprior;         % offset point numbers by already existing point indices.
  [M,ret]=addobject(Ko,M,'curve',no,0); 
end;
end;

fclose(fid);
return;

%-----------------------------------------------------------
function M=set_selectables(Ko,M,co,data)
switch co,
  case 'points',
%      disp('point');
      val = get(Ko.typefields(1),'userdata');
      if val==0, 
	 switch data(1),
	  case 0,
	  M.selectablepoints = find(M.pointtypes==0);  Kmessage('all Points');
	  case 1,
	  M.selectablepoints = find(M.pointtypes==1 & M.pointreferences==data(2));  
	  Kmessage('Free Vertices only');
	  case 2,   
	  M.selectablepoints = find(M.pointtypes==2);  Kmessage('Curve Vertices only');
	  case 3,
	  M.selectablepoints = find(M.pointtypes==3);  % surface vertices
	 end;
       else 
	 M.selectablepoints=[];
       end;
       set(Ko.typefields(1),'userdata',~val);
       M = point_indexing(Ko,M,0);
       
  case 'curves',  
      val = get(Ko.typefields(2),'userdata');
      val=~val;
      set(Ko.typefields(2),'userdata',val);
      M = curve_indexing(Ko,M,0);
      
  case 'surfaces',  
      Kmessage('not implemented');
 end;
    
 for k=1:2,
   v = get(Ko.typefields(k),'userdata');
   if v==0,
     set(Ko.typefields(k),'backgroundcolor','w');
   else
     set(Ko.typefields(k),'backgroundcolor','y');
   end;
 end;
 return;

%-----------------------------------------------------------
function M=point_indexing(Ko,M,co,num)
  switch co,
   case 'up',    
    if isempty(M.selectablepoints), Kmessage('no points to select'); return; end;
    M = select_point(Ko,M,0);
    npts = length(M.selectablepoints);
    M.selectedptindex = min(M.selectedptindex+1,npts);
    M.redpoint=M.selectablepoints(M.selectedptindex); 
    M = select_point(Ko,M,M.redpoint);    
    return;
   case 'down',
    if isempty(M.selectablepoints), Kmessage('no points to select'); return; end;
    M = select_point(Ko,M,0);
    npts = length(M.selectablepoints),
    M.selectedptindex = max(M.selectedptindex-1,1);
    M.redpoint=M.selectablepoints(M.selectedptindex); 
    M = select_point(Ko,M,M.redpoint);    
   case 'set',
    M.lastpoint=M.redpoint;
    M = select_point(Ko,M,0);
    M.redpoint=num;
    Kmessage(sprintf('red: %i  last: %i  [redvertx: %i lastvertx %i]',...
		     M.redpoint,M.lastpoint,...
		     find(M.vertextable==M.redpoint),...
		     find(M.vertextable==M.lastpoint)));
    
    M = select_point(Ko,M,M.redpoint);
   case 'splinept',
    M.lastpoint=M.redpoint;
    M = select_point(Ko,M,0);
    M.redpoint=num;
    M = select_point(Ko,M,M.redpoint);
    M.splinepointindex = M.splinepointindex+1;
    M = memorize_obj(Ko,M,'point',M.splinepointindex);    % put point in cell. 
  end;
set(Ko.selects(1),'string',sprintf('%i',M.redpoint));
return;

%-----------------------------------------------------------
function M=curve_indexing(Ko,M,co,num)
switch co,
   case 'up',
     M = select_point(Ko,M,0);
     M = select_curve(M,0);
     ncur=length(M.curvehandles);
     M.redcurve = min(ncur,M.redcurve+1);
     M = select_curve(M,M.redcurve);
   case 'down',
     M = select_point(Ko,M,0);
     M = select_curve(M,0);
     ncur=length(M.curvehandles);
     M.redcurve = max(0,M.redcurve-1);
     M = select_curve(M,M.redcurve);   
   case 'set',
     M = select_point(Ko,M,0);
     M = select_curve(M,0);
     M.lastcurve=M.redcurve;
     M.redcurve = num;
     M = select_curve(M,M.redcurve);
 end;
 set(Ko.selects(2),'string',sprintf('%i',M.redcurve)); 
 return;
 
%-----------------------------------------------------------
function M=select_point(Ko,M,indx)
% if a point is selected, it's relevant parameters must be in the
% numerical display (xyz-fields). So we need to program the 
% xyz editor appropriately.
if indx==0,
  if M.redpoint ~= 0 ,
     h = M.pointhandles(M.redpoint);
%     disp([h,M.redpoint,ishandle(h)]);
     [msize,colo,mark]=objproperties('point',M.pointtypes(M.redpoint),0);
     set(h,'color',colo,'markersize',msize,'marker',mark);
     set(h,'hittest','on');  % attente: that could create a problem if points should not be selectable
     M.redpoint=0;
   end;
else,
   h = M.pointhandles(indx);
   [msize,colo,mark]=objproperties('point',M.pointtypes(indx),1);
   set(h,'color',colo,'markersize',msize,'marker',mark);
   set(h,'hittest','off');
   M.redpoint = indx;
   M = display_numerics(Ko,M);
 end;
set(Ko.selects(1),'string',sprintf('%i',M.redpoint)); 
return;

%-----------------------------------------------------------
function M= display_numerics(Ko,M)
% display the coordinates in the xyz editor.
red = M.redpoint;
if red==0; return; end;
typus = M.pointtypes(red);
a = M.xyz(red,:);
setedit_values(Ko.xyzedits,a);
return;

%----------------------------------------------------------- 
function M=select_curve(M,indx)
 if indx==0,
   if M.redcurve,
      h = M.curvehandles(M.redcurve);
      [msize,colo,mark]=objproperties('curve',M.curvetypes(M.redcurve),0);
      set(h,'color',colo,'linewidth',msize);
    end;
 else,
    h = M.curvehandles(indx);
    [msize,colo,mark]=objproperties('curve',M.curvetypes(indx),1);
    set(h,'color',colo,'linewidth',msize);
    M.redcurve=indx;
 end;
  
  if M.redcurve > 0, 
    M.selectablepoints = M.curves{M.redcurve};
    M.selectedptindex=0;
  else,
    M.selectablepoints=[];
    M.selectedptindex=0;
  end;
   
return;
 
%-----------------------------------------------------------
function [msize,colo,mark]=objproperties(objecttype,subtype,onoff)
if onoff == 0,
  switch objecttype,
    case 'point',
    switch subtype,
      case 0,  msize=6; colo = 'b'; mark='.';   % data point
      case 1,  msize=18; colo= 'g'; mark='.';   % free vertex
      case 2,  msize=18; colo= 'm'; mark='.';   % curve point vertex
      case 3,  msize=16; colo= 'w'; mark='.';  % Bezier curve control vertex on an edge
      case 4,  msize=6; colo= 'w'; mark='diamond';   % Bezier surface control vertex.
    end;
    case 'curve',
    switch subtype,
      case 0,  msize=0.5;  colo='b'; mark='none'; % data curve;
      case 1,  msize=1;  colo='c'; mark='none'; % constructed curve
    end;
  end;
elseif onoff==1,	% selected mode.  
  switch objecttype,
    case 'point',
    switch subtype,
      case 0,  msize=20; colo = 'r'; mark='.';   % data point
      case 1,  msize=28; colo= 'r'; mark='.';   % free vertex
      case 2,  msize=8; colo= 'y'; mark='o';   % curve point
      case 3,  msize=10; colo= 'y'; mark='hexagra';  % Bezier curve control vertex on an edge
      case 4,  msize=10; colo= 'y'; mark='diamond';   % Bezier surface control vertex.
    end;
    case 'curve',
    switch subtype,
      case 0,  msize=1;  colo='r'; mark='none'; % data curve;
      case 1,  msize=2;  colo='y'; mark='none'; % constructed curve
    end;
  end;
end;

%-----------------------------------------------------------  
function M=turn_into_vertex(Ko,M)
red = M.redpoint;
if red>0,
  M = select_point(Ko,M,0);
  M.pointtypes(red)=1;
  h = M.pointhandles(red);
  xyz=[get(h,'xdata'),get(h,'ydata'),get(h,'zdata')];
  [msize,colo,mark]=objproperties('point',1,0); 
  [M,nindex]= addobject(Ko,M,'point',xyz,1);
  M = select_point(Ko,M,nindex);
else,
  Kmessage('no red point defined');
end;
return;
%-----------------------------------------------------------

function M = upgrade_edges_table(Ko,M)
% change booo done
% this function goes through the list of edges and finds
% the ones that are only defined as straightedges. If there 
% is now a curve in the drawing that connects the endpoints of such
% an edge, it will replace the straightedge. 
[nc,n3]=size(M.edgetable);
disp('upgrade_edges');
for k=1:nc,
  if M.edgetable(k,1)==0,             % if not registered as a curve already.
    nv=abs(sort(compact(M.A,k)));
    nv1=nv(1);
    nv2=nv(2);
    np1=M.vertextable(nv1);
    np2=M.vertextable(nv2);      
    % next, go throught all curves and find one that connects np1 and np2.
    for l=1:length(M.curves);
      ncu = M.curves{l};
      n1 = find(ncu==np1);
      n2 = find(ncu==np2);
      if n1*n2>0, 
	M.edgetable(k,1:3)=[l,np1,np2];         % make proper entry [curve, point1 point2]
	delete(M.edgehandles(k));               % erase straight line edge,
	M.edgehandles(k)=0;                     % clear handle of the straight edge
	disp([nv1,nv2,np1,np2,n1,n2]);
      end;
    end;
  end;
end;
%-----------------------------------------------------------

function M=popout_cube(Ko,M)
%  gets the memorized points and pops out a cube.
  mem=get_memory_content(M);
  if isempty(mem),
    return;
  end;
  if length(mem)<8,
    mem = [mem,zeros(1,8-length(mem))];
  end;
  
  if ~all(mem(1:4)) 
    Kmessage('Points 1-4 not defined ');
    return;
  end;
  
  no = mem;
  x14 = M.xyz(mem(1:4)',1:3);
  idef = find(mem(5:8));
  ndef = length(idef); % number of defined ones.
  udef = find(mem(5:8)==0);  % undefined.
  
  if ndef ~= 0,
    defs = mem(5:8);
    defs = defs(idef);   % defined points
    x58=x14;
    x58(idef,1:3) = M.xyz(defs,1:3);
    dx = x58-x14;
    dx = sum(dx)/ndef;
    if length(udef)>0,
      x58(udef,1) = x58(udef,1) + dx(1);
      x58(udef,2) = x58(udef,2) + dx(2);
      x58(udef,3) = x58(udef,3) + dx(3);
    end;
  else,
    x58 = x14;
    dx = cross(x14(2,1:3)-x14(1,1:3),x14(4,1:3)-x14(1,1:3))+...
         cross(x14(3,1:3)-x14(2,1:3),x14(1,1:3)-x14(2,1:3))+...
         cross(x14(4,1:3)-x14(3,1:3),x14(2,1:3)-x14(3,1:3))+...
         cross(x14(1,1:3)-x14(4,1:3),x14(3,1:3)-x14(4,1:3));
    dx = dx/4;
    nn=norm(dx);
    dx = dx/sqrt(nn);
    x58(:,1) = x58(:,1)+dx(1);
    x58(:,2) = x58(:,2)+dx(2);
    x58(:,3) = x58(:,3)+dx(3);
  end;
  
  for k=1:length(udef),
    [M,np]= addobject(Ko,M,'point',x58(udef(k),1:3),1);   % is going to be a vertex.
    M.vertextable = [M.vertextable,np];
    no(udef(k)+4)= np;
  end;
  for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
  name = get(Ko.namefield,'string')
  if isempty(name) | isblank(name), name='xxxxxxxx'; end;
  name = make_a_label(name,8);
  M = add_a_cell(Ko,M,'cube',no,name);
  M=clear_object_memory(Ko,M,'points');
  return;
%-----------------------------------------------------------

function M=popout_prism(Ko,M)
%  gets the memorized points and pops out a cube over a triangle
  mem=get_memory_content(M);
  if isempty(mem),
    return;
  end;
  if length(mem)<6,             % append zeros for empty points
    mem = [mem,zeros(1,6-length(mem))];
  end;
  
  if ~all(mem(1:3))             % min 3 points needed
    Kmessage('Points 1-3 not defined ');
    return;
  end;
  
  no = mem;
  x13 = M.xyz(mem(1:3)',1:3);
  idef = find(mem(4:6));
  ndef = length(idef); % number of defined ones.
  udef = find(mem(4:6)==0);  % undefined.
  
  if ndef ~= 0,
    defs = mem(4:6);
    defs = defs(idef);   % defined points
    x36=x13;
    x36(idef,1:3) = M.xyz(defs,1:3);
    dx = x36-x13;
    dx = sum(dx)/ndef;
    if length(udef)>0,
      x36(udef,1) = x36(udef,1) + dx(1);
      x36(udef,2) = x36(udef,2) + dx(2);
      x36(udef,3) = x36(udef,3) + dx(3);
    end;
  else,
    x36 = x13;
    dx = cross(x13(2,1:3)-x13(1,1:3),x13(3,1:3)-x13(1,1:3))+...
         cross(x13(3,1:3)-x13(2,1:3),x13(1,1:3)-x13(2,1:3))+...
         cross(x13(1,1:3)-x13(3,1:3),x13(2,1:3)-x13(3,1:3));
    dx = dx/3;
    nn=norm(dx);
    dx = dx/sqrt(nn);
    x36(:,1) = x36(:,1)+dx(1);
    x36(:,2) = x36(:,2)+dx(2);
    x36(:,3) = x36(:,3)+dx(3);
  end;
  
  for k=1:length(udef),
    [M,np]= addobject(Ko,M,'point',x36(udef(k),1:3),1);   % is going to be a vertex.
    M.vertextable = [M.vertextable,np];
    no(udef(k)+3)= np;
  end;
  for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
  name = get(Ko.namefield,'string')
  if isempty(name) | isblank(name), name='xxxxxxxx'; end;
  name = make_a_label(name,8);
  M = add_a_cell(Ko,M,'prism',no([4 5 2 1 6 3]),name);
  M=clear_object_memory(Ko,M,'points');
  return;
%-----------------------------------------------------------

function set_tissue_namefield(Ko,tissuename)
  % used by the program maketissuemenu thru the case 'tissuename'
    
  set(Ko.namefield,'string',tissuename);
  return;
  
function M=popout_roof(Ko,M)
%  gets the memorized points and pops out a roof prism
  mem=get_memory_content(M);
  if isempty(mem),
    return;
  end;
  if length(mem)<6,
    mem = [mem,zeros(1,6-length(mem))];
  end;
  
  if ~all(mem(1:4)) 
    Kmessage('Points 1-4 not defined ');
    return;
  end;
  
  no = mem;
  x14 = M.xyz(mem(1:4)',1:3);
  idef = find(mem(5:6));
  ndef = length(idef); % number of defined ones.
  udef = find(mem(5:6)==0);  % undefined.
  
  if ndef ~= 0,
    defs = mem(5:6);
    defs = defs(idef);   % defined points
    x56=[0.5*(x14(1,1:3)+x14(2,1:3));
         0.5*(x14(3,1:3)+x14(4,1:3))];
    x56(idef,1:3) = M.xyz(defs,1:3);   % overwrite if defined.
%    dx = x56-x14;
%    dx = sum(dx)/ndef;
    dx = x56(1,1:3) - 0.5*(x14(1,1:3)+x14(2,1:3));
    if length(udef)>0,
      x56(udef,1) = x56(udef,1) + dx(1);
      x56(udef,2) = x56(udef,2) + dx(2);
      x56(udef,3) = x56(udef,3) + dx(3);
    end;
  else,
    x56=[0.5*(x14(1,1:3)+x14(2,1:3));      % start with interpolations
         0.5*(x14(3,1:3)+x14(4,1:3))];     % between 1-2 and 3-4
    % calculate vector shift perpendicular to plane 1 2 3 4
    dx = cross(x14(2,1:3)-x14(1,1:3),x14(4,1:3)-x14(1,1:3))+...
         cross(x14(3,1:3)-x14(2,1:3),x14(1,1:3)-x14(2,1:3))+...
         cross(x14(4,1:3)-x14(3,1:3),x14(2,1:3)-x14(3,1:3))+...
         cross(x14(1,1:3)-x14(4,1:3),x14(3,1:3)-x14(4,1:3));
    dx = dx/4;
    nn=norm(dx);
    dx = dx/sqrt(nn);
    x56(:,1) = x56(:,1)+dx(1);
    x56(:,2) = x56(:,2)+dx(2);
    x56(:,3) = x56(:,3)+dx(3);
  end;
  
  % if the points x56 are indeed generated above, they need to be put in the vertex table:  
  for k=1:length(udef),
    [M,np]= addobject(Ko,M,'point',x56(udef(k),1:3),1);   % is going to be a vertex.
    M.vertextable = [M.vertextable,np];
    no(udef(k)+4)= np;
  end;
  for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
  name = get(Ko.namefield,'string')
  if isempty(name) | isblank(name), name='xxxxxxxx'; end;
  name = make_a_label(name,8);
  M = add_a_cell(Ko,M,'prism',no,name);
  M=clear_object_memory(Ko,M,'points');
  return;
%-----------------------------------------------------------


function M=start_object(Ko,M)
    M=clear_object_memory(Ko,M,'points');    
    M=begin_spline(Ko,M);  % understand that also other things might be generated that way..
return;

%-----------------------------------------------------------

function M=create_object(Ko,M,subc)
switch subc,
  case 'spline',
  % get the points from memory: 
  mem = get_memory_content(M)  % fetch the buttons
  if isempty(mem),
    M=begin_spline(Ko,M);    % what about some feedback here: turn red the button
    set(findobj('tag','Splinestarter'),'backgroundcolor','r');
    return;
  end;
  [M,flg]= addobject(Ko,M,'curve',mem(:),1,0);   % the equivalent of end_spline.
  M = upgrade_edges_table(Ko,M);
  set(findobj('tag','Splinestarter'),'backgroundcolor',[0.7 0.7 0.7]);
  % added later: 
  if M.splinepointindex>0,  % clear memory
    M=clear_object_memory(Ko,M,'points');
  end;
end;
return;
%-----------------------------------------------------------
function M=make_curve_point(Ko,M)
redc = M.redcurve;
if redc ~= 0,
  M = select_point(Ko,M,0);
  xyz=get_curvepoints(M,redc);
  Y = cubicsplinecurve(xyz,1,0);     % generate beginning point.
  Y = [Y(1),Y(2),Y(3),0];
  [M,nindex]=addobject(Ko,M,'point',Y,2,redc);  % where does it get it's curve from?
   M = select_point(Ko,M,nindex);
else,
  Kmessage('No red curve to make curve point');
end;   
%-----------------------------------------------------------  
function Kmessage(str,rec)
if nargin==1,
  rec=0;
end;
o = findobj('tag','FEMKonstructorMessage');
set(o,'string',str);
%x=get(o,'backgroundcolor');
%set(o,'backgroundcolor','y'); drawnow;
%set(o,'backgroundcolor',x);   drawnow;
us = get(o,'userdata');
if ~isempty(us),
  if (us==1 & rec==1) | us==2,
    fid=fopen('Kmessagefile.rec','a');
    fprintf(fid,'%s\n',str);
    fclose(fid);
  end;
end;

return;
%-----------------------------------------------------------  
function Kmessage_recorder(onoff)
%  A callback to a radiobutton , see construct.m  
  o = findobj('tag','FEMKonstructorMessage');
  switch onoff,
   case 1
    set(o,'userdata',1);
   case 2
    set(o,'userdata',2);
   case 0,
    set(o,'userdata',[]);    
    fid=fopen('Kmessagefile.rec','a');
    fprintf(fid,'------------------------------------\n');
    fclose(fid);
  end
return;
%-----------------------------------------------------------  
function xyz=get_curvepoints(M,ncurve)
np = M.curves{ncurve};
for k=1:length(np),
  h = M.pointhandles(np(k));
  xyz(k,:) = [get(h,'xdata'),get(h,'ydata'),get(h,'zdata')];
end;

%-----------------------------------------------------------
function M = toggle_autoproject(Ko,M,onoff)
if onoff==0,
  M.autoupdate=1;
elseif onoff == 1,
  M.autoupdate=2;
end;
%-----------------------------------------------------------
function M = move_to_projection(Ko,M,method)
% GET point from ebenen master and move selected point there
% method - 0: project onto plane, 1 move parallel to plane.
if nargin<3, method=0; end;

red=M.redpoint;
if red==0, return; end;
ret = ebenenmaster('clicked');  % see ebenenenmaster.m for the returned data structure
check=ret{1};    % tells if plane was hit
xyz= ret{2};     % projected point
pars=ret{3};     % parametric description of intersection
nor = ret{4};    % nor normal direction of the plane
if check==1,
  if method==1,
    xa = M.xyz(red,1:3) ;                 % old value
    xyz =( [(xa-xyz')*nor]*nor+ xyz)';    % leave normal component unchanged.
  end;
  M=change_xyz_values(Ko,M,xyz);
  if M.autoupdate>0,
    M= update_hierarchy(Ko,M);
  end;
end;
return;

%-----------------------------------------------------------
function xyz = get_projection_coords(Ko,M);
ret = ebenenmaster('clicked');  % see ebenenenmaster.m for the returned data structure
check=ret{1};    % tells if plane was hit
xyz= ret{2};     % projected point
pars=ret{3};     % parametric description of intersection
nor = ret{4}; 
if check==0,
  xyz=[];
end;
return;

%-----------------------------------------------------------
 
function M = project_curve_inplane(Ko,M)
redc = M.redcurve;
if redc==0, return; end;
cu = M.curves{redc};
for k=1:length(cu),
  npt = cu(k);
  typus = M.pointtypes(npt);
  if typus>0,
    h = M.pointhandles(npt);
    xyz=M.xyz(npt,1:3);
    xyz = ebenenmaster('project',xyz);  % ask the guy to project onto currently active plane
    M.xyz(npt,1) = xyz(1);
    M.xyz(npt,2) = xyz(2);
    M.xyz(npt,3) = xyz(3);
    M=change_xyz_values(Ko,M,xyz,npt);
  end;
end;
M = update_hierarchy(Ko,M);
return;

%-----------------------------------------------------------

function M = project_point_inplane(Ko,M)
% If a point is selected it is projected into the currently active plane (see
% ebenenmaster).
red=M.redpoint;
if red==0, return; end;
xyz = M.xyz(red,1:3);
xyz = ebenenmaster('project',xyz);      % ask the guy to project onto currently active plane
M=change_xyz_values(Ko,M,xyz);
if M.autoupdate>0,
  M= update_hierarchy(Ko,M);
end;

return;

%-----------------------------------------------------------

function M = change_xyz_values(Ko,M,xyz,red)
% function called if a point's coordinates are modified in
% the GUI or by other means. 
 if nargin == 2, 
   xyz=getedit_values(Ko.xyzedits);
 end;
 if nargin < 4,
   red = M.redpoint;
   if red == 0,   return; end;
 end;
 typus = M.pointtypes(red);
 h = M.pointhandles(red);
 if M.autoupdate==2,
   xyz(1:3) = ebenenmaster('project',xyz(1:3)');  % allways project onto current plane.
 end;
switch typus
  case 0,               % normal point. 
   Kmessage('Data point selected - can''t be changed');
  case 1,               % vertex 
   M.xyz(red,1:3) = [xyz(1),xyz(2),xyz(3)];
   set(h,'xdata',xyz(1),'ydata',xyz(2),'zdata',xyz(3));
   necke = find(M.vertextable==red);
   rerender_bezEcke(Ko,M,necke); 

  case 2,              % curve vertex (who's coordinates are changed by the r paramter)
   nc = M.pointreferences(red);    % what curve is it on
   npts = M.curves{nc};            % point indices of curve
   Xcontr = M.xyz(npts(:),1:3);    % coordinates of those points.
   x = cubicsplinecurve(Xcontr,1,xyz(4)); % take r field as parameter.
   M.xyz(red,1:3) = x(:)';
   M.xyz(red,4) = xyz(4);
   set(h,'xdata',x(1),'ydata',x(2),'zdata',x(3));
   setedit_values(Ko.xyzedits(1:3),x(1:3));
   
  case 3,     % edge control vertex.
   nedge = M.pointreferences(red);  % number of edge this point is control vertex of.
   M.xyz(red,1:3) = [xyz(1),xyz(2),xyz(3)];
   set(h,'xdata',xyz(1),'ydata',xyz(2),'zdata',xyz(3));
%   Kmessage(sprintf('edgepoint of edge %i',nedge));
   %% redo that edge and attached surfaces.
   M=rerender_bezedge(Ko,M,nedge);

  case 4,   
   nsrf=M.pointreferences(red);  
   M.xyz(red,1:3) = [xyz(1),xyz(2),xyz(3)];
   set(h,'xdata',xyz(1),'ydata',xyz(2),'zdata',xyz(3));
%   Kmessage(sprintf('point of surface %i',nsrf));
   %% redo that surface
   rerender_bezier(Ko,M,nsrf);
end;
return;

%-----------------------------------------------------------
function M= memorize_obj(Ko,M,subcode,data)
switch subcode,
  case 'point',
    nmem = data;
    if M.redpoint==0,
      Kmessage('no point to memorize');
      return;
    end;
    M.memory(nmem) = M.redpoint;
    set(Ko.merker(nmem),'backgroundcolor','y');  
  case 'curve',
    nmem = data;
    if M.redcurve==0,
      Kmessage('no curve to memorize');
      return;
    end;
    
    M.curvememory(nmem) = M.redcurve;
    set(Ko.curvemerker(nmem),'backgroundcolor','y');  
    
    case 'loop',
    Kmessage('not implemented');
  end;

%-----------------------------------------------------------  
function indx = recall_points(Ko,M)
indx = M.memory(find(M.memory));
return;
%-----------------------------------------------------------
function M = clear_object_memory(Ko,M,type)
switch type,
case 'points',
 for k=1:16,
  %  M.memory(k)=0;
  set(Ko.merker(k),'backgroundcolor','w');  
 end;
 M.memory=[];
 M.memoryperm=[];
 M.splinepointindex=0;
 
case 'curves',
 for k=1:12,
   set(Ko.curvemerker(k),'backgroundcolor','w');
 end;
 M.curvememory=[];
 case 'loops',
  for k=1:12,
   set(Ko.loopmerker(k),'backgroundcolor','w');
 end;
 M.loopmemory=[];
 disp('not implemented');
end;
return;
%-----------------------------------------------------------
function M = show_memory(Ko,M,onoff)
% label with a little number.
if ~isempty(M.memoryperm),
  Mmem = M.memory(M.memoryperm);
else,
  Mmem = M.memory;
end;
memlen = length(find(Mmem));

if onoff==1,
  for k=1:memlen,
    str = sprintf('%i',k);
    npt = Mmem(k);
    xyz = M.xyz(npt,1:3);
    h = text(xyz(1), xyz(2), xyz(3)+1, str);
    set(h,'color','g');
    pth = M.pointhandles(npt);
    co = get(pth,'color');
    set(pth,'visible','off');  % (color','w');
    fu.color= co;
    fu.number = npt;  % index
    set(h,'userdata',fu,'tag','Tempory_Number');
  end;
elseif onoff==0,
  fl = findobj('tag','Tempory_Number');
  if isempty(fl), return; end;
  for k=1:length(fl),
    fu = get(fl(k),'userdata');
    npt=fu.number;
    pth=M.pointhandles(npt);
    set(pth,'visible','on', 'color',fu.color);  % restore color
    delete(fl(k));
  end;
end;
%-----------------------------------------------------------
function M = permute_memory(Ko,M)    % uses edit field tagged editmemorysequence
 ed = findobj('tag','editmemorysequence');
 str = get(ed,'string');
 no = length(M.memory);      % that's how many integers I will try to read from the edit/
 [perm,nx]=sscanf(str,'%i',no);
 if nx == no,
   M.memoryperm = perm';
 elseif nx==0,
   M.memoryperm =[];
 else
   Kmessage('Wrong number of permutation indices');
 end;
return
%-----------------------------------------------------------
function Mmem=get_memory_content(M)
% delivers a set of point indices.
if ~isempty(M.memoryperm),
  Mmem = M.memory(M.memoryperm);
else,
  Mmem = M.memory;
end;
return;
%-----------------------------------------------------------
function M=zoominout(Ko,M,val)
red = M.redpoint;
dxa=10;
if val==0,
  ax = Ko.axisrange;
  if ~isempty(ax), axis(ax); end;
else,
  if red==0, 
    Kmessage('No red point selected to zoom to');
  else,
    xyz = M.xyz(red,1:3);
    Ko.axisrange=axis;
    ax = [xyz(1)-dxa, xyz(1)+dxa, xyz(2)-dxa, xyz(2)+dxa, xyz(3)-dxa, xyz(3)+dxa];
    axis(ax);
  end;
end;
%-----------------------------------------------------------
function M=save_picture(Ko,M,filename)
fid = fopen(filename,'w');
M = reorder_objects(M);
for k=1:length(M.pointtypes),
  h = M.pointhandles(k);
  typ = M.pointtypes(k);
  ref = M.pointreferences(k);
  xyz = M.xyz(k,:);
  fprintf(fid,'P %2i %5i   %5.2f %5.2f %5.2f    %5.2f %5.2f %5.2f\n',typ,ref,xyz);
end;

for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  np = length(no);
  fprintf(fid,'C %2i %5i %4i  ', typ, ref, np);
  for j=1:np,
    fprintf(fid,'%4i ',no(j));
  end;
  fprintf(fid,'\n');
end;
fprintf(fid,'E\n');
fclose(fid);
return;

%-----------------------------------------------------------
function M=save_macro(Ko,M,filename)
% This function assumes that the whole design is a macro and saves it correspondingly to a
% file. The file can later be included to generate a macro cell that may be deformed and
% moved before being glued together with the rest of the system. 
% 

nameregistration('append',filename);

fid = fopen(filename,'w');
M = reorder_objects(M);
for k=1:length(M.pointtypes),
  h = M.pointhandles(k);
  typ = M.pointtypes(k);
  ref = M.pointreferences(k);
  xyz = M.xyz(k,:);
  switch typ,
    case 0, let ='D';
    case 1, let ='V';
    case 2, let ='D';
    case 3, let ='I';
    case 4, let ='D';
  end;
  fprintf(fid,'%c %3i  %5.2f %5.2f %5.2f\n',let,k,xyz(1:3));
end;

% write out the curves as well. 
for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  np = length(no);
  fprintf(fid,'S  %4i', np);
  for j=1:np,
    fprintf(fid,'%4i ',no(j));
  end;
  fprintf(fid,'\n');
end;
% write out the cells.
for k=1:length(M.cells),
  typ = M.cells{k}{1};
  no  = M.cells{k}{2};
  switch typ, 
    case 'cube', o='C';
    case 'prism', o='P';
    case 'tetra', o='T';
  end;
  fprintf(fid,'%s ',o);
  for j=1:length(no), fprintf(fid,'%4i ',no(j)); end; 
  fprintf(fid,'\n');
end;

for k=1:length(M.cells),
  clab=M.celllabels{k};
  if isequal(clab,'        '),
    clab= make_a_label(filename,8);
  end;
  fprintf(fid,'L %4i %8s\n',k,clab);
end;

fprintf(fid,'E\n');
fclose(fid);
return;

%-----------------------------------------------------------
function M=save_selected_pict(Ko,M,filename)
% looks for the splines that have a type greater than 0
% and writes only the points out that are referenced by the curves.
% 
tagfield = zeros(length(M.pointhandles),1);
for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  np = length(no);
  if typ>0,
    for j=1:np,
      tagfield(no(j))=1;
    end;
  end;
end;
% give each referenced point a new index:
ta = find(tagfield==1);  
nta = length(ta);
tagfield(ta) = [1:nta];

fid = fopen(filename,'w');
% write out only the points that are used in the 
% curves of type > 0:
for k=1:length(M.pointtypes),
  h = M.pointhandles(k);
  typ = M.pointtypes(k);
  ref = M.pointreferences(k);
  if tagfield(k)~=0,
    typ=1; ref=1;
    xyz = M.xyz(k,:);
    xyz(4:6) = [0 0 0];
    fprintf(fid,'P %2i %5i   %5.2f %5.2f %5.2f    %5.2f %5.2f %5.2f\n',typ,ref,xyz);
  end;
end;
% write out only the curves with type>0 (while modifying the point indices)

for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  if typ>0,
    ref=M.curvereferences(k);
    ref=1;
    np = length(no);
    no = tagfield(no);  % replace indices.
    fprintf(fid,'C %2i %5i %4i  ', typ, ref, np);
    for j=1:np,
      fprintf(fid,'%4i ',no(j));
    end;
    fprintf(fid,'\n');
  end;
end;
fprintf(fid,'E\n');
fclose(fid);
return;
%-----------------------------------------------------------

function [M,flg]=load_picture(Ko,filename,flag)
if flag==0,
  M = makeGeometry;
end;
fid = fopen(filename,'r');
if fid<0, 
  flg=0;
  Kmessage(['Could not open ',filename(1:30)','...']);
  return;
end;
eof=0;
while ~eof,
  [C,n] = fscanf(fid,'%c',1);
  switch C,
  case 'E',
  eof=1;
  break;
  case 'P',
%  [no,n]=fscanf(fid,'%i',1);
  [kenn,n]=fscanf(fid,'%i',2);
  typ = kenn(1);
  ref = kenn(2);
  [xyz,n]=fscanf(fid,'%f',6);
  [M,ret]= addobject(Ko,M,'point',xyz,typ,ref);
%  disp(xyz');

   case 'C',
  [kenn,n] = fscanf(fid,'%i',2);
  typ = kenn(1);
  ref = kenn(2);
  [np,n] = fscanf(fid,'%i',1);  % number of points
  [no,n] = fscanf(fid,'%i',np); % indices
%  disp(no');
  [M,ret]=addobject(Ko,M,'curve',no,typ,ref); 
end;
end;
% after this, redefine curve points.
% new: Pressing prism/tetra/cube will generate a loose object
%of that type, given current point selections in memory. 
%Pressing clr will just discard it again. 
%Pressing take will implement it: put it in the general topology. 
%Putting it in the general topology, means: modifying A,B,C and
%generating constructed curves between the points as well as
%storing of the control vertices. Hence we need the following
%new data structures%
%
%macrovertices:  an array of point indices that are used as the 
%                vertices in the topological structure
%macrocurves     an array of curve indices that are used as the
%                edges in the topological structure
%
%		etc. 
%              
%Later we'll probably need more crap.
% 
%-----------------------------------------------------------
function [M,flg]= show_topology(Ko,M,subcode)
% turn memorypoints into tetra, prism or cube
% subcode is 'tetra','cube',or 'prism'

% first look it there is already an object like that:

flg=1;
killit=0;

switch subcode,
case 'tetra',     np=4;
case 'prism',     np=6;
case 'cube',      np=8;
case 'none',      killit=1;
case 'take',      [M,flg] = make_a_topology(Ko,M);  return;  % Assumes you got already a tempory one.
end;

h = findobj('tag','temporary_topological_object');
if ~isempty(h),
  if ishandle(h), 
    delete(h);
  end;
end;
if killit, return; end;
mem = get_memory_content(M);  % fetch the buttons

if length(mem) ~= np,
  Kmessage(['Inappropriate number of points in memory for ',subcode]);
  flg=0;
  return;
end;

[a,b,c] = Glue([],[],[],subcode,0,[1:np]);

[ma,na] = size(a);
xyz = zeros(ma*2+(ma-1),3);    % reserver for edges
kk=0;
for k=1:ma,
  ll = compact(a,k);
  ll = abs(ll);
  m = mem(ll);
  kk=kk+1; xyz(kk,:) = M.xyz(m(1),1:3);  % get first point's coods
  kk=kk+1; xyz(kk,:) = M.xyz(m(2),1:3);  % get 2nd point's coords
  if k<ma, kk=kk+1; xyz(kk,:) = [NaN,NaN,NaN]; end; % empty filler (don't draw);
end;
h = line('xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3),'color','y');
userdat.type = subcode; 
userdat.indx = mem;

set(h,'tag','temporary_topological_object','userdata',userdat);
% this way the information is stored in the tempory object. make_a_topology uses this information.
return;
 
%----------------------------------------------------------- 
function [M,flg]= make_a_topology(Ko,M)
% makes a topological model permanent.
flg=1;
%
% first look if there is a tempory one:
tmp = findobj('tag','temporary_topological_object');
if isempty(tmp),
  flg=0; 
  Kmessage('Can''t find a topological element');
  return;
end;

% So we got one: find out what it is; then kill it.
usd = get(tmp,'userdata'); 
mem = usd.indx  ;             % the points of which it's made

delete(tmp);
Kmessage(['Ok - I''ll make a real ',usd.type]);

% What does it take:
% A the vertices
% First I need to get indices for the vertices in question: 
% For that we ask another program to  give me the vertex indices.
% If one or more of these points are not in the vertex list,
% They are going to be added to it.
% 
[M,vindex]=vertex_indices(M,mem);    
% 
% remember: vindex contains indices into M.vertextable. In that table we
% also have only indices which then point to the actual data (M.xyz). 
[nb,mb]=size(M.B);
[M.A,M.B,M.C,newedges]=Glue(M.A,M.B,M.C, usd.type, 0, vindex);  % extern
[nbn,mbn]=size(M.B);
newloops = [(nb+1):1:nbn];    % (loops are added at the end of the matrix).
% So we just got the new topology. Now we still need to create the edges.
% newedges  contains the row indices that where added to the A matrix.
% So we create a different representation of that:
tab = compact(M.A,newedges);  % a cell array
% this contains the edges that need to be drawn. Since they have not been
% in the topology before (Glue returns only new lines), they all need to 
% be plotted, from the initial point to the end. 
[M,flg]=make_edges(Ko,M,tab);
% this first of all just draws curves, but gives them a special index. 
[M,flg]=make_loops(Ko,M,newloops);
return;

%Each loop or surface is associated with additional information.  There may be additional
%points in the surface, and there may be additional handles for lines or other methods to
%draw the surface. For now I will keep an array of handles, one handle for each surface.
% If the surface has no additional internal points then it has to be drawn without them:
% this is done by linear interpolation of the splines. 

%% end make topology
%-----------------------------------------------------------
function [M,flg]=make_edges(Ko,M,tab)
% tab is a cell containing signed vertex indices.
flg=1;
if ~iscell(tab), 
  tab={tab};
end;
for k=1:length(tab),
  ll = tab{k};
  n1 = abs(min(ll));    % beginning point index stripped sign
  n2 = max(ll);         % end point index.
% n1 and n2 are indices in M.vertextable
% np1 = M.vertextable(n1);
% np2 = M.vertextable(n2);
% np1 and np2 are indices in M.xyz
% This has to change. Replace by a Kurvenstueck if it exists. 
  [M,ncurve] = edges_table(Ko,M,n1,n2);  
%  [M,ncurve] = addobject(Ko,M,'curve',[np1,np2],1);   % curve of type 1. 
  if ncurve==0, Kmessage('make_edges got it all wrong: see addobject'); flg=0; end;
end;
  
%-----------------------------------------------------------

function [M,vdx]=vertex_indices(M,inx)
% inx: a bunch of point indicies that need to be a vertex
nhighest = length(M.vertextable);
for k=1:length(inx),
  iv = find(~(M.vertextable - inx(k)));  % is in table?
  if isempty(iv),             
     nhighest=nhighest+1;
     M.vertextable(nhighest) = inx(k); % if not, put in table. 
     iv = nhighest;
  end;
  vdx(k) = iv;
end;
return;
%-----------------------------------------------------------
function [M,flg]=make_loops(Ko,M,newloops)
% newloops are the indices (into B) of the new loops. 
flg=1;
for lo=newloops,
 nl = compact(M.B,lo);
 M.loopintern{lo}{1}=nl;
 M.loopintern{lo}{2}=[];   % place holder for point references. 
end;
%-----------------------------------------------------------
function M = update_hierarchy(Ko,M,allways)
% called to update the curves after points have changed. 
% 
% the redpoint may have changed, and all curves that refere to it,
% must change too.
if nargin==2,
  allways=0;
  if (M.redpoint==0 & M.redcurve==0), return; end;  % OK, nothing at all to do
end;

for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  f = find(no==M.redpoint);
  if ~isempty(f) | allways==1,
    h = M.curvehandles(k);
    xyz = M.xyz(no,1:3);
    [n,n3]=size(xyz);
    xyz = cubicsplinecurve(xyz,n*5,'equal');
    set(h,'xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3));
    % now check if a curve point is defined relative to this curve:
    M=update_curvepoint(Ko,M,k);
  end;
end;
% finally , update all straightlines that might have changed:
[ma,na]=size(M.A);
for k=1:ma,
  h = M.edgehandles(k);
  if h~=0 & ishandle(h),
    nv = compact(M.A,k);
    nv1 = abs(min(nv));
    nv2 = abs(max(nv));
    nn=M.vertextable([nv1,nv2]);
    X = [M.xyz(nn(1),1:3);M.xyz(nn(2),1:3)];
    set(h,'xdata',X(:,1),'ydata',X(:,2),'zdata',X(:,3));
  end;
end;
return;

%-----------------------------------------------------------
function M = update_dependence(Ko,M,npoint)
% called to update the curves after a  
% point has changed. 
% the redpoint may have changed, and all curves that refere to it,
% must change too.

for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  f = find(no==npoint);
  if ~isempty(f),
    h = M.curvehandles(k);
    xyz = M.xyz(no,1:3);
    [n,n3]=size(xyz);
    xyz = cubicsplinecurve(xyz,n*5,'equal');
    set(h,'xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3));
    % now check if a curve point is defined relative to this curve:
    M=update_curvepoint(Ko,M,k);    % attention this may result in recursion.
  end;
end;
return;

%-----------------------------------------------------------
function M=update_curvepoint(Ko,M,kcurve);
ncpt = find(M.pointtypes==2 & M.pointreferences==kcurve);
if ~isempty(ncpt),
  for k=1:length(ncpt),
    red = ncpt(k);
    xyz = M.xyz(red,:);
    h = M.pointhandles(red);
    npts = M.curves{kcurve};                  % point indices of curve
    Xcontr = M.xyz(npts(:),1:3);    % coordinates of those points.
    x = cubicsplinecurve(Xcontr,1,xyz(4)); % take r field as parameter.
    M.xyz(red,1:3) = x(:)';
    M.xyz(red,4) = xyz(4);
    set(h,'xdata',x(1),'ydata',x(2),'zdata',x(3));
    M = update_dependence(Ko,M,red);
  end;
end;
return;

%-----------------------------------------------------------
function [M,flg]= insert_curve_point(Ko,M)
% On the currently selected spline, use the currently selected 
% moving point (curve point) and make it one of the spline's 
% control points. 
flg=0;
if isempty(M.redpoint) | isempty(M.redcurve) 
  Kmessage('no curve or no curve point selected'); 
  return; end;
if M.redcurve == 0, Kmessage('no curve selected'); return; end;
if M.redpoint == 0, Kmessage('no point selected'); return; end;
if M.pointtypes(M.redpoint) ~= 2, Kmessage('no curve point selected'); return; end;

% OK jetzt ham wir Ruhe: Kurve und Punkt selektiert:
pts = M.curves{M.redcurve};

xyz = M.xyz(M.redpoint,:);
r = xyz(4);                 % parameter value.
if r==0, 
  newpts = pts; Kmessage('Selected curve point too close to start of spline'); return;
elseif r==1,
  newpts = pts; Kmessage('Selected curve point too close to end of spline'); return;
else
  Y = M.xyz(pts,:);
  u = pathparameters(Y);      % parameterized way through points. 
  uu = [r;u(:)];
  [u,ix]=sort(uu);
  newpts = [M.redpoint,pts(:)'];  % 
  newpts = newpts(ix);        % reorder according to parameter value.
end;
 
redp = M.redpoint;             % remember redpoint
M = select_point(Ko,M,0);     % unselect it
M.pointtypes(redp) = 1;       % 
M.curves{M.redcurve} = newpts;
M = select_point(Ko,M,redp);  % select it again.
M = update_hierarchy(Ko,M); 
flg=1;
return;

%--------------------------------------------------------
function M=reorder_objects(M);
[ptyp ,O ]= sort(M.pointtypes);
M.pointhandles = M.pointhandles(O);
M.xyz = M.xyz(O,:);
M.pointtypes=ptyp;
M.pointreferences = M.pointreferences(O);
[io,inperm]=sort(O);   % inverse permutaion
for k=1:length(M.curves),
  pts = M.curves{k};
  pts = inperm(pts);
  M.curves{k} = pts;
end;
if M.redpoint~=0,
  M.redpoint = inperm(M.redpoint);
end;
return;
%-----------------------------------------------------------
function M=show_spline_direction(Ko,M,val)
% shows a little line segment to recognise beginning of curve'
bu = findobj('tag','splinedirection');
if ~isempty(bu) & ishandle(bu),
  delete(bu);
end;

if val==1,
  if M.redcurve ~=0, 
    nc=M.redcurve;
    pts = M.curves{nc};
    Xcontr = M.xyz(pts(:),1:3);                 % coordinates of those points.
    x = cubicsplinecurve(Xcontr,11,[0.0:0.01:0.1]); % take r field as parameter.
    bu = line('xdata',x(:,1),'ydata',x(:,2),'zdata',x(:,3),...
	'linewidth',5,'color','r');
    set(bu,'tag','splinedirection');
  end;
end;
%-----------------------------------------------------------

function [M,flg] = attach_pt_to_spline(Ko,M,how);   
% subcode 1: at end, subcode 0: at beginning
% verify first if a spline is selected and if
% there is at least one point in the memory.
% 
  flg=0;
  if M.redcurve == 0, Kmessage('no curve selected'); return; end;
  
  mem = get_memory_content(M)
  
  if isempty(mem) | mem==0, Kmessage('no points in memory'); return; end;
  
  pts =M.curves{M.redcurve};
  
  if how==0, pts = [mem(:); pts(:)]';
  elseif how==1, pts = [pts(:); mem(:)]';
  end;
  
  M.curves{M.redcurve} = pts;
  M = upgrade_edges_table(Ko,M);
  M = update_hierarchy(Ko,M,1); 
  M=clear_object_memory(Ko,M,'points'); 
  flg=1;
  return;
  
%
%-----------------------------------------------------------

function M = join_connected_splines(Ko,M)
  % two consecutively selected splines are jointed, the previously
  % selected one is attached to the currently selected one.
  % In any case, redcurve remains the selected curve.

  if M.redcurve==0
    Kmessage('No curve selected');
    return;
  end;
  if M.lastcurve==0,
    Kmessage('Only one curve selected');
    return;
  end;
  % Now assume there are two selected:
  pts1 = M.curves{M.redcurve};  n1 = length(pts1);
  pts2 = M.curves{M.lastcurve}; n2 = length(pts2);   % last curve is
                                                     % previous curve

% four possibilities of joining:
  pts=[];
  
  if pts1(n1)==pts2(1),
    pts = [pts1;pts2(2:n2)];
  elseif pts1(n1) == pts2(n2),
    pts = [pts1; flipud(pts2(1:n2-1))];
  elseif pts1(1)==pts2(1),
    pts = [flipud(pts1); pts2(2:n2)];
  elseif pts1(1)==pts2(n2),
    pts = [flipud(pts1); flipud(pts2(1:n2-1))];
  end;

  disp(pts);
  
  if isempty(pts),
    Kmessage('two selected curves have no common end point');
    return;
  end;
  
  kilc=M.lastcurve;
  M.lastcurve=0;
  redmem=M.redcurve;
  M=destroy_current_curve(Ko,M,kilc);   % att: will reset redcurve
  M.redcurve=redmem;
  M.curves{M.redcurve} = pts;
%  M = update_hierarchy(Ko,M); 
  M = upgrade_edges_table(Ko,M);
  M = update_hierarchy(Ko,M,1); 
  return;
%-----------------------------------------------------------
  
function M=set_preselect(Ko,M,subcode)
% subcode: either 'points' or 'curves';
pts = findobj('tag','point');
cur = findobj('tag','curve');

M=surface_indexing(Ko,M,'off');   % an afterthought, see the disfunct Bezier stuff

switch subcode,
  case 'points',
  for k=1:length(cur), set(cur(k),'hittest','off'); end;
  for k=1:length(pts), 
    uu = get(pts(k),'userdata');
    if uu(2)>=1, 
      set(pts(k),'hittest','on','buttondownfcn',...
	  'FEMKonstructor(''hit'',''point'',get(gco,''userdata''));','color','y','visible','on'); 
    end;
  end;
  case 'curves',
  for k=1:length(pts), set(pts(k),'hittest','off'); end;
  for k=1:length(cur),
    set(cur(k),'hittest','on','buttondownfcn',...
	'FEMKonstructor(''hit'',''curve'',get(gco,''userdata''));'); 
  end;
end;
return;

%-----------------------------------------------------------

function M=begin_spline(Ko,M)
% in variation of the above subroutine we define a callback for all points.
  pts = findobj('tag','point');
  cur = findobj('tag','curve');
  for k=1:length(cur), set(cur(k),'hittest','off'); end;
  for k=1:length(pts), 
    uu = get(pts(k),'userdata');
    if uu(2)>=1, 
      set(pts(k),'hittest','on','buttondownfcn','FEMKonstructor(''hit'',''splinept'',get(gco,''userdata''));'); 
    end;
  end;
  M.splinepointindex = 0;
  return;
  
%-----------------------------------------------------------
function M = hit_object(Ko,M,subcode,data)
% data is the userdata, extracted in backup
%
switch(subcode),
case 'point', 
 M = point_indexing(Ko,M,'set',data(1));
case 'splinept', 
 M = point_indexing(Ko,M,'splinept',data(1));
case 'curve', 
 M = curve_indexing(Ko,M,'set',data(1));
case 'surface',
 M = surface_indexing(Ko,M,'pick',data); 
end;
return;
%-----------------------------------------------------------

function [M,nedg]=edges_table(Ko,M,nv1,nv2)
% MODIFY: MAKE BEZIER SPLINE PIECES BY CALCULATING THREE INTERMEDIATE POINTS
% 
% finds the first  Kurvenstuck (segment of a spline) that connects n1 and n2.
% if there isn't any, it will make a straight line. 
% M.edgetable is 0 if there is a straight line, and contains the
% curve index if the edge is a curve segment. 
% The information of the curve segment is retrieved by 
% looking up the indices of the vertices of the segment in the A matrix
% then looking up the point indices corresponding to the vertices in
% M.vertextable and finally by finding these indices in the curve
ncurve=0;
%fprintf(1,'make edge between %i and %i\n',nv1,nv2);

np1 = M.vertextable(nv1);    % corresponding point numbers. 
np2 = M.vertextable(nv2);    % corresponding point numbers. 
nedg=0;

for k=1:length(M.curves),
  ncu = M.curves{k};
  n1 = find(ncu==np1);
  n2 = find(ncu==np2);
  if n1*n2 > 0,                % if yes, curve connectes the points
    [ncurve,n3]=size(M.edgetable);
    ncurve=ncurve+1;
    M.edgetable(ncurve,1:3) =[k,np1,np2];     % now k is curve number
    [ncurve,n3]=size(M.edgetable);
    M.edgehandles(ncurve) = 0;
    nedg=ncurve;      % to flag that found a curve piece.
    return;
  end;
end;
% otherwise if no connecting curve could be found make a straight line:
[nedg,n3]=size(M.edgetable);
nedg=nedg+1;
M.edgetable(nedg,1:3) = [0,np1,np2];
[nedg,n3]=size(M.edgetable);
X = [M.xyz(np1,1:3);M.xyz(np2,1:3)];
h = line('xdata',X(:,1),'ydata',X(:,2),'zdata',X(:,3),'color','g');
set(h,'tag','straightedge');
set(h,'hittest','off');         % initial 
M.edgehandles(nedg) = h;
return;

%-----------------------------------------------------------
function [M,flg]= split_curve_at_pt(Ko,M);   % 
% grabs curve from memory, and uses redpoint to split
% the curve if the redpoint is on the curve.
flg=0;
np=M.redpoint;
if isempty(np) | np==0,
  Kmessage('no point selected');
  return;
end;
nc = M.curvememory;
if isempty(nc) | length(nc)>1,
  Kmessage('none or more than one curve selected');
  return;
end;

np=np(1);
nc=nc(1);
c = M.curves{nc};
pp = find(c==np);
if isempty(pp), 
  Kmessage('Split point not part of spline');
  return;
end;
flg=1;

nc1 = c(1:pp);
nc2 = c(pp:end);

typ = M.curvetypes(nc);
ref = M.curvereferences(nc);
M.curves{nc}=nc1;
[M,flg]= addobject(Ko,M,'curve',nc2,typ,ref);
M = update_hierarchy(Ko,M); 
return;
%-----------------------------------------------------------

% To proceed: I put a cube on the screen and move the vertices around, 
% I put another on the screen and glue it by selecting correspondences and so on.
% The glueing is done by selecting corresponding points and clicking a button. 
% 
%-----------------------------------------------------------

function [M,flg]=create_free_cell(Ko,M,type,name)
flg=1;
Kmessage(sprintf('Create a free %s at current settings',type));
switch type,
  case 'cube', npoints=8; 
        xyz=[0 0 0; 1 0 0; 1 1 0; 0 1 0;0 0 1; 1 0 1; 1 1 1; 0 1 1];
  case 'prism', npoints=6;
        xyz=[0 0 0; 1 0 0; 1 1 0; 0 1 0;0.5 0 1; 0.5 1 1];
  case 'tetra', npoints=4;
        xyz=[0 0 0; 1 0 0; 1 1 0; 0.5 0.5 1];
end;

% set editor for geometry to decent innitial values:
% keep old position but clear rotation and resize
parms = getedit_values(Ko.topoedits);
parms(4:9) = [0 0 0 10 10 10];
setedit_values(Ko.topoedits,parms);
% 
for k=1:npoints,
 [M,np]=addobject(Ko,M,'point',10*xyz(k,:),1);
 M.vertextable = [M.vertextable,np];
 vindex(k) = length(M.vertextable);
end;
[nb,mb]=size(M.B);
[M.A,M.B,M.C,newedges]=Glue(M.A,M.B,M.C, type, 0, vindex);  % extern
[nbn,mbn]=size(M.B);
newloops = [(nb+1):1:nbn];    
% (loops are added at the end of the matrix).
% So we just got the new topology. Now we still need to create the edges.
% newedges  contains the row indices that where added to the A matrix.
% So we create a different representation of that:
tab = compact(M.A,newedges);  % a cell array
% this contains the edges that need to be drawn. Since they have not been
% in the topology before (Glue returns only new lines), they all need to 
% be plotted, from the initial point to the end. 
[M,flg]=make_edges(Ko,M,tab);
% this first of all just draws curves, but gives them a special index. 
[M,flg]=make_loops(Ko,M,newloops);

ncell=length(M.cells);
M.cells{ncell+1} = {type, vindex};
M.celllabels{ncell+1} = name;
M.freecell=ncell+1;
M.freecellnumber=1;
[M,flg]=move_free_cell(Ko,M);
return;
%-----------------------------------------------------------
function [M,flg]=move_free_cell(Ko,M)
% The free cell, for now, should only be the last cell generated. 
% (easy extension: the last N cells to make macros possible).
% Changes:
% This program was modified from one that can only move one cell
% at a time, namely the last cell, to one that can move and deform
% multiple cells. 
% 
ncell = M.freecell;   % first free cell.
flg=0;
if ncell==0, flg=0; Kmessage('no free cell'); return; end;


if M.freecellnumber==1 & M.specialcell==0,
  type = M.cells{ncell}{1};
  nvertx = M.vertextable(M.cells{ncell}{2});   % vertex numbers-> point numbers.
else,
  type='multi';  % for below's switch   
  nvertx = [];   % find all vertices associated with the macro.:
% just put all point references that can be found in a long array and
% then sort out multiple references. 
  for k=1:M.freecellnumber,
    nv = M.vertextable(M.cells{ncell+k-1}{2});
    nvertx = [nvertx;nv(:)];
  end;
  nmaxspli = M.macrosplines(1);    % first macro spline
  nsplines = M.macrosplines(2);    % nunber of splines in macro.
  for k=1:nsplines,
    cu = M.curves{nmaxspli+k};
    nvertx = [nvertx;cu(:)];
  end;
  mx = max(nvertx);    % remove multiple references....
  aa = zeros(mx,1);
  aa(nvertx)=1;
  nvertx = find(aa);   % list of verticies in macro each vertex only once.
end;
    
parms = getedit_values(Ko.topoedits);
xm = parms(1); ym=parms(2); zm=parms(3); 
rx = parms(4); ry=parms(5); rz=parms(6);
a = parms(7); b = parms(8); c=parms(9);

xmid = [xm,ym,zm];
sca = [a,b,c];

sx = sin(rx*pi/180);
cx = cos(rx*pi/180);
Rx = [ 1   0    0;
       0  cx  -sx;
       0  sx   cx];

sz = sin(rz*pi/180);
cz = cos(rz*pi/180);
Rz = [ cz   sz   0;
       -sz   cz   0;
       0    0    1];

sy = sin(ry*pi/180);
cy = cos(ry*pi/180);
Ry = [cy  0  -sy;
      0   1    0;
      sy  0   cy];

RR = Rz*Ry*Rx;  % rotation matrix.

% nps= M.vertextable(nvertx);     % actual points
hs = M.pointhandles(nvertx);
nps=nvertx;

switch type,
case 'cube',
  npts=8;
  xyz(1,:) = [-a,-b,-c]/2;
  xyz(2,:) = [+a,-b,-c]/2;
  xyz(3,:) = [+a,+b,-c]/2;
  xyz(4,:) = [-a,+b,-c]/2;
  xyz(5,:) = [-a,-b,+c]/2;
  xyz(6,:) = [+a,-b,+c]/2;
  xyz(7,:) = [+a,+b,+c]/2;
  xyz(8,:) = [-a,+b,+c]/2;

  case 'prism',
  npts=6;
  xyz(1,:) = [-a,-b,-c]/2;
  xyz(2,:) = [+a,-b,-c]/2;
  xyz(3,:) = [+a,+b,-c]/2;
  xyz(4,:) = [-a,+b,-c]/2;
  xyz(5,:) = [0, -b, c]/2;
  xyz(6,:) = [0, +b, c]/2;
  
  case 'tetra',
  npts=4;
  xyz(1,:) = [-a,-b,-c]/2;
  xyz(2,:) = [+a,-b,-c]/2;
  xyz(3,:) = [0 ,+b,-c]/2;
  xyz(4,:) = [0,0,c]/2;
  
  case 'multi',
  if isempty(M.macrocoords),     % get initial state 
    xyz = M.xyz(nps,1:3);
    M.macrocoords=xyz;
    newm=1;
  else,
    xyz=M.macrocoords;
    newm=0;
  end
  [npts,n3]=size(xyz);
  xmin = min(xyz);
  xmax = max(xyz);
  xmid0=(xmin+xmax)/2;
  siz=(xmax-xmin)/2;
  if newm~=0,
    vals = [xmid0(:)',[0 0 0],siz(:)'];
    setedit_values(Ko.topoedits,vals); 
  else
    xyz = xyz-ones(npts,1)*xmid0; % move point of gravity to orinin
    xyz = xyz.*(ones(npts,1)*1*(sca./siz));  
  end;
  
  otherwise,
  disp('no cell of that type');
  return
end;   % switch ende

xyz = (RR*xyz')';
xyz = xyz+ones(npts,1)*[xm,ym,zm];

for k=1:length(hs),
  M.xyz(nps(k),1:3) = [xyz(k,1),xyz(k,2),xyz(k,3)];
  set(hs(k),'xdata',xyz(k,1),'ydata',xyz(k,2),'zdata',xyz(k,3)); %%  ,'color','y');
end;
M = update_gluons(M);
M = update_hierarchy(Ko,M,1);
return;

%-----------------------------------------------------------
function M = show_cell_numbers(Ko,M,onoff,mactype)
  if nargin<4, mactype=[]; end;
  if nargin<3, onoff=1; end;
  o = findobj('tag','cellcentertext');
  if ~isempty(o), delete(o); end;
  if isempty(onoff), onoff=1; end;
  
  if onoff==0, return; end;
  switch onoff,
   case 1,
    for k=1:length(M.cells);
      nif=1;
      if ~isempty(mactype),
	nif=isequal(mactype,M.celllabels{k});
      end;
      if nif==1,
	no = M.cells{k}{2};
	no = M.vertextable(no);
	labl=M.celllabels{k};
	dr = tissueMaster('draw?',labl,k);
	if dr==1,
	  xm = mean(M.xyz(no,1:3));   % point of gravity
	  h = text(xm(1),xm(2),xm(3),sprintf('%i',k));
	  set(h,'color','m','fontsize',18,'tag','cellcentertext');
	end;
      end;
    end;
   case 2,
    k = mactype;
    no = M.cells{k}{2};
    no = M.vertextable(no);
    xi = M.xyz(no,1:3);
    xm = mean(xi);   % point of gravity
    h = text(xm(1),xm(2),xm(3),sprintf('%i',k));
    set(h,'color','m','fontsize',18,'tag','cellcentertext');
    for j=1:length(no),
      h = text(xi(j,1),xi(j,2),xi(j,3),sprintf('%i',j));
      set(h,'color','r','fontsize',20,'tag','cellcentertext');
    end;
  end;
  return
  
%-----------------------------------------------------------      
function M = show_vertices_numbers(Ko,M,onoff)
% label with a little number.

nvtx = M.vertextable;
hs = M.pointhandles(nvtx);
xyz = M.xyz(nvtx);
npts = length(nvtx);

if onoff==1,
  for k=1:npts,
    pth = hs(k);
    vis = get(pth,'visible');
    if isequal(vis,'on'),
      str = sprintf('%i',k);
      xyz = M.xyz(nvtx(k),1:3);
      h = text(xyz(1), xyz(2), xyz(3), str);
      set(h,'color','g','fontsize',18);
      co = get(pth,'color');
      set(pth,'visible','off');  % (color','w');
      fu.color= co;
      fu.number = nvtx(k);  % point index
      set(h,'userdata',fu,'tag','Tempory_Number');
    end;
  end;
  
elseif onoff==0,
  fl = findobj('tag','Tempory_Number');
  if isempty(fl), return; end;
  for k=1:length(fl),
    fu = get(fl(k),'userdata');
    npt=fu.number;
    pth=M.pointhandles(npt);
    set(pth,'visible','on', 'color',fu.color);  % restore color
    delete(fl(k));
  end;
end;

%-----------------------------------------------------------
function [M,flg]=glue_near_blocks(Ko,M)
% Take all vertices and check if any are very close together. 
% if one can find two that are close together, call once glue_point_to_point.
flg=0;
threshold=1;
[check,n1,n2]=find_close_points(M,threshold);

while check,
  [M,flg]=glue_point_to_point(Ko,M,n1,n2);
  [check,n1,n2]=find_close_points(M,threshold);
end;
return;
%%%%%%%%%%%%%%%%%%
function [check,n1,n2]=find_close_points(M,threshold)
% returns vertex indices in n1 and n2.
nvert = length(M.vertextable);
check=0;
for k=1:nvert-1,
  np1 = M.vertextable(k);
  n1=k;
  for j= k+1:nvert;
    np2 = M.vertextable(j);
    n2 = j;
    dif = diff(M.xyz([np1,np2]',1:3));
    dis = sum(dif.^2);
    if dis<threshold,
      check=1; return;
      break;
    end;
  end;
end;
return;

%-----------------------------------------------------------
function M = newpoint_between_two(Ko,M,ninter)
% takes the last selected and the currently selected points and includes the
% edge that connects them, a new point. 
if nargin==2, ninter=1; end;

if M.redpoint==0 | M.lastpoint==0, flg= -1; return; end;
n2 = find(M.vertextable==M.redpoint);
if isempty(n2), 
  Kmessage('Last point selected for glueing operation is not vertex - no action');
  flg=-1;
  return;
end;

n1 = find(M.vertextable==M.lastpoint);
if isempty(n1), 
  Kmessage('Previous point selected for glueing operation is not vertex - no action');
  flg=-1;
  return;
end;

%disp([n1,n2]);

np1 = M.vertextable(n1);
np2 = M.vertextable(n2);

% make intermediate points.
x1 = M.xyz(np1,1:3);
x2 = M.xyz(np2,1:3);

xnew = sum(M.xyz([np1,np2],1:3))/2;   % intermediate point coordinates
if M.autoupdate==2, 
  xnew = ebenenmaster('project',xnew(:)');  % automatically project on active plane if needed. 
end;

%  find the curve that connects them: Condition they have to be neighbors:
for kcur=1:length(M.curves),
  cu = M.curves{kcur}; 
  cu=cu(:);
  i1=find(cu==np1);
  i2=find(cu==np2);
  if ~isempty(i1) & ~isempty(i2) & abs(i1-i2)==1,    % si ils sont voisins,
    dt = 1/(ninter+1);
    nindex=zeros(ninter,1);   
    xip = curvepiece(M.xyz,cu,np1,np2,ninter+2);    % clculate n+2 equally spaced points, take n
    xip = xip(2:ninter+1,1:3);                      % take the n middle ones only
    for k=1:ninter,
      xi = xip(k,:);
      if M.autoupdate==2, 
	xi = ebenenmaster('project',xi(:)');  % automatically project on active plane if needed. 
      end;
      [M,nindex(k)]=addobject(Ko,M,'point',xi,3,kcur); %% generate intermediate points.
    end;
    len=length(cu);
%   disp([len,i1,i2]);
    if i1<i2, 
      nucu = [cu(1:i1)', nindex(:)',  cu(i2:len)'];
    elseif i2<i1,
      nindex=flipud(nindex);
      nucu = [cu(1:i2)', nindex(:)', cu(i1:len)'];
    end;
    M.curves{kcur}=nucu;
    M=update_hierarchy(Ko,M);
    break;  % don't look further if points are created.
  end;
end;
return;


%-----------------------------------------------------------  
function  [M,flg]= glue_point_to_point(Ko,M,n1,n2)
% what does it do: two points where highlighted one after the other,
% join them into one if they are vertices. 
% 
flg=1;

if nargin==2,
  if M.redpoint==0 | M.lastpoint==0, flg= -1; return; end;
  n2 = find(M.vertextable==M.redpoint);
  if isempty(n2), 
    Kmessage('Last point selected for glueing operation is not vertex - no action');
    flg=-1;
    return;
  end;
  
  n1 = find(M.vertextable==M.lastpoint);
  if isempty(n1), 
    Kmessage('Previous point selected for glueing operation is not vertex - no action');
    flg=-1;
    return;
  end;
end;
% So, according to the above n1 and n2 are vertex indices.
% The second one selected moves to the first one selected, c'est a dire
% redpoint moves to lastpoint, c'est a dire
% n2 moves to n1.
if isempty(n1) | isempty(n2),
  Kmessage('Points selected for glueing operation are not vertex - no action');
  return;
end;

if M.gluons == 1,
  M=make_gluon(M,n1,n2);    % postpone glueing and create a "gluon": a line that connects the two
  return;
end;

if ~isempty(n1) & ~isempty(n2),
 [M.A,M.B,M.C,rmpts,rmedges,rmloops] = Klebzusammen(M.A,M.B,M.C,[n1,n2]);  % glue-together
end;

h = M.edgehandles(rmedges);
ta = M.edgetable(rmedges,1);
for k=1:length(h),
  if (ta(k)==0),   delete(h(k)); end;  
end;
M.edgehandles(rmedges) =[];
M.edgetable(rmedges,:)=[];
% 
[M,flg]=replace_points(Ko,M,n2,n1);
M = update_hierarchy(Ko,M,1);
return;

%-----------------------------------------------------------
function M=make_gluon(M,n1,n2)
% This function is called after it has been established that n1 and n2 are vertices
% It generates a thin line between the two points and hangs the information into it that 
% is needed to later glue the two points together. Since the indexing can changed, we have 
% to rely on information stored in the graphic handles. 

np1 = M.vertextable(n1);
np2 = M.vertextable(n2);
h1 = M.pointhandles(np1);
h2 = M.pointhandles(np2);
xyz = M.xyz([np1,np2]',1:3);
hgluon = line('xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3),...
    'tag','gluon','linestyle','-.','color','y','hittest','off');
set(hgluon,'userdata',[h1,h2]);
return;
%-----------------------------------------------------------

function show_gluons(M,how,col)
gl = findobj('tag','gluon');
if isempty(gl), return; end; 

switch how,
case 'color',
  for k=1:length(gl),
    has=get(gl(k),'userdata');
    h1=has(1);
    h2=has(2);
    set(h1,'color',col);
    set(h2,'color',col);
  end;
case 'invisible',
set(gl,'visible','off');
case 'visible',
set(gl,'visible','on');
end;
return;

%-----------------------------------

function M=update_gluons(M)
% update the position of the gluons if changes of xyz coordinates have taken place.
% 
gl = findobj('tag','gluon');
if isempty(gl), return; end; 

for k=1:length(gl),
  has=get(gl(k),'userdata');
  h1=has(1);
  h2=has(2);
  n1 = get(h1,'userdata'); n1 = n1(1);
  n2 = get(h2,'userdata'); n2 = n2(1);
  xyz = M.xyz([n1,n2]',1:3);
  set(gl(k),'xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3));
end;
return;

%---------------------------------------------------------------------------
function M = glue_gluon_pairs(Ko,M)
%
glu = findobj('tag','gluon');   % find '-.' lines and represent a glue rule.
for k=1:length(glu),
  has = get(glu(k),'userdata');
  h1=has(1);
  h2=has(2);
  n1 = get(h1,'userdata'); n1 = n1(1);
  n2 = get(h2,'userdata'); n2 = n2(1);
  n1 = find(M.vertextable==n1);
  n2 = find(M.vertextable==n2);
  [M.A,M.B,M.C,rmpts,rmedges,rmloops] = Klebzusammen(M.A,M.B,M.C,[n1,n2]);
  h = M.edgehandles(rmedges);
  ta = M.edgetable(rmedges,1);
  for j=1:length(h),
    if (ta(j)==0),   delete(h(j)); end;  
  end;
  M.edgehandles(rmedges) =[];
  M.edgetable(rmedges,:)=[];
  % 
  [M,flg]=replace_points(Ko,M,n2,n1);
  M = update_hierarchy(Ko,M,1);
  delete(glu(k));
end;
return;



%-----------------------------------------------------------
%%% end of 

%%
% State: Need to fix the problem with the edges list: 
% Should not contain point numbers. 
%% 
%----------------------------------------------------------------------
function [M,flg]=replace_points(Ko,M,n2,n1)
% n1 and n2 are vertex indices not point indicies.
% This is a function that needs some work: Free points are just killed, 
% points that are part of a spline need to update the spline. 
% if a point is part of more than one spline, it can not be killed...
% n2 moves to n1. If n2 is on a spline (or multiples) 
% it means it has to take the coordinates
% of n1 and the spline(s) have to be updated. However, if
% n2 is on a straight edge, it means the straight edge has to be updated.
% 
% where there used to be n1, there is now n2, that is n2 copies the
% information from n1 while n1 get's distroyed. The information associated with 
% n1 is it's location
%n2 is a point in a spline, let's say the 3rd one. now this becomes n1, that
%is at the third position of the spline that contains n2, n1 is stored.
%
%Now let's say n2 is point in an edge: then it means that where n2 is stored,
%n1 will be stored. Is that simple or what?
%Then point n2 will be destoyed since it's no longer refereing to anything
% (I wrote this stupid comment since I got all confused for some time. 
% first check multiple refereed points, remove them from the kill list.
np2 = M.vertextable(n2);       % points' indices
np1 = M.vertextable(n1);       % replace by point indices.
flg=0;
%forall splines that contain n2 , replace n2 index by n1 index. end;
%forall edges that contain n2 replace n2 by n1. 
%forall cells that conatin n2 replace n2 by n1 done in reorder_points below.
for k=1:length(M.curves),
  cu = M.curves{k};
  nse = find(cu==np2);
  if ~isempty(nse), 
    cu(nse)=np1;
    M.curves{k} = cu;
  end;
end;

[ne,n3]=size(M.edgetable);
for k=1:ne,
  if M.edgetable(k,2) == np2, M.edgetable(k,2) = np1; end;
  if M.edgetable(k,3) == np2, M.edgetable(k,3) = np1; end;
end;

% now that all references to n2 have been erased, we can erase n2:
M = select_point(Ko,M,0);
h = M.pointhandles(np2);
if ishandle(h), delete(h);end;
M.pointhandles(np2)=0;
M = reorder_points(M,n2,n1);     
% cells also need to be updated via the pthand thingi.
M = update_hierarchy(Ko,M,1);
return;
%-------------------------------------------------

function M=kill_point(Ko,M,nkill,update)
% can kill a control point of a spline but not a really vertex. 
% 
% first look if it's in the vertext list:
if nargin==2,
  nkill=M.redpoint;
  if nkill==0; return; end;
  M = select_point(Ko,M,0);   % unselect it.
end;
if nargin<4, update=1; end;

ii = find(M.vertextable==nkill); 
if ~isempty(ii),
  Kmessage('This point cannot be removed: would change topology');
  return;
end;
h = M.pointhandles(nkill);
if ishandle(h), delete(h); end;
M.pointhandles(nkill)=0;
pthand = 1:length(M.pointhandles);
pthand(M.pointhandles==0) = 0;
fo = find(pthand);
nleft = length(fo);
g = 1:nleft;
pthand(fo)=g;
for k=1:length(M.curves),
  cu = M.curves{k};
  cu = pthand(cu);
  M.curves{k}=cu(find(cu));
end;
nonul = find(M.pointhandles);
M.pointhandles=M.pointhandles(nonul);
M.pointreferences=M.pointreferences(nonul);
M.pointtypes = M.pointtypes(nonul);
M.xyz = M.xyz(nonul,:);
for k=1:length(M.pointhandles),
  nus=get(M.pointhandles(k),'userdata');
  nus(1) = k;
  set(M.pointhandles(k),'userdata',nus);
end;
if update==1,
  M = update_hierarchy(Ko,M,1);
end;
return;


%-------------------------------------------------


function M=reorder_points(M,n2,n1)
% this function detects empty slots in the point list, and reorders all.
% first create a mapping array:
pthand = 1:length(M.pointhandles);
pthand(M.pointhandles==0) = 0;


fo = find(pthand);
nleft = length(fo);
g = 1:nleft;
pthand(fo)=g;
for k=1:length(M.curves),
  cu = M.curves{k};
  cu = pthand(cu);
  M.curves{k}=cu;
end;
mo=M.vertextable;

oldnumvertex = length(M.vertextable);

M.vertextable = pthand(M.vertextable);  % exclude cleared vertex -> zero
M.vertextable = M.vertextable(find(M.vertextable));  % eliminate zeros

% reorder vertex references in the cells. tricky tricky!!.
% fprintf('n1=%i  n2=%i',n1,n2);

pt = zeros(oldnumvertex,1);
if n1>n2;
pt(n2)=n1-1;
else
pt(n2)=n1;
end;

[pts,remap]=reordering_array(pt);

for k=1:length(M.cells),
  no = M.cells{k}{2};
  M.cells{k}{2}=remap(no);
end;

M.edgetable(:,2) = pthand(M.edgetable(:,2));
M.edgetable(:,3) = pthand(M.edgetable(:,3));

%e1=M.edgetable;

%disp([e0,e1]);


nonul = find(M.pointhandles);
M.pointhandles=M.pointhandles(nonul);
M.pointreferences=M.pointreferences(nonul);
M.pointtypes = M.pointtypes(nonul);
M.xyz = M.xyz(nonul,:);
for k=1:length(M.pointhandles),
  nus=get(M.pointhandles(k),'userdata');
  nus(1) = k;
  set(M.pointhandles(k),'userdata',nus);
end;
return;
%-----------------------------------------------------------

function M=change_curve_color(Ko,M,col)
if M.redcurve==0, Kmessage('no curve selected, no color change'); return; end;
h = M.curvehandles(M.redcurve);
M = select_curve(M,0);
M.redcurve=0;
set(h,'color',col);
return;


function M=destroy_current_curve(Ko,M,kilc)
% If kilc is not defined, 
% the currently selected curve is killed and references to it are replaced.
% otherwise it is assumed that kilc is the index of the curve to be destroyed.
  
  if nargin<3, 
    kilc=M.redcurve;
    M.redcurve=0;
  end;

  if kilc==0, return; end;

  if M.curvehandles(kilc) & ishandle(M.curvehandles(kilc)), 
    delete(M.curvehandles(kilc)); 
  end;
  nc = length(M.curvehandles);
  M.curvehandles(kilc) =[];
  M.curvetypes(kilc)=[];
  M.curvereferences(kilc)=[];
  
  map = ones(nc,1);
  map(kilc)=0;
  map = cumsum(map);
  map(kilc)=0;
  M.curves{kilc}={};
  kk=0;
  for k=1:length(M.curves),
    if ~isempty(M.curves{k}), 
      kk=kk+1;
      ncu{kk} = M.curves{k};
    end;
  end;
  M.curves = ncu;
  fo=find(M.edgetable(:,1));
  M.edgetable(fo,1) = map(M.edgetable(fo,1));   
  M = rerender_graph(Ko,M);
  M = upgrade_edges_table(Ko,M);
  return;

%-----------------------------------------------------------
function [pt,map]=reordering_array(cond)
% create intermediate address mapping fields for 
% reordering indicies after the delection of certains
% cond contains for each non-zero entry by which index
% the current index should be replaced. 
% For example, cond = [0 1 0 0 3 0 0 0 0 4]' generates:
% 
%    cond   pt   map
%     0     1     1
%     1     0     1
%     0     2     2
%     0     3     3
%     3     0     3
%     0     4     4
%     0     5     5
%     0     6     6
%     0     7     7
%     4     0     4
% Compris?
% 

cond = cond(:);
pt = [1:length(cond)]';
pt(find(cond)) = 0;
fo = find(pt);
nleft = length(fo);
g = [1:nleft]';
pt(fo)=g;
map=pt;
map(find(cond~=0)) = cond(find(cond~=0));
%disp([cond,pt,map]);
return;

%-----------------------------------------------------------
function M = rerender_graph(Ko,M)
% called to rerender the whole thing after its either just changed or 
% reloaded from file. 
[npt,nx]=size(M.xyz);

for k=1:npt,  
  typ=M.pointtypes(k);
  ref=M.pointreferences(k);
  [msize,colo,mark]=objproperties('point',typ,0);
  if M.pointhandles(k) & ishandle(M.pointhandles(k)),
    set(M.pointhandles(k),'xdata',M.xyz(k,1), 'ydata',M.xyz(k,2), 'zdata',M.xyz(k,3));
    set(M.pointhandles(k),'userdata',[k,typ]);
  else
    M.pointhandles(k) = line('xdata',M.xyz(k,1), 'ydata',M.xyz(k,2), 'zdata',M.xyz(k,3));
  end;
  set(M.pointhandles(k),'marker',mark,'markersize',msize,'color',colo,'tag','point','userdata',[k,typ]);
end;

for k=1:length(M.curves),
  no = M.curves{k};
  typ=M.curvetypes(k);
  ref=M.curvereferences(k);
  h = M.curvehandles(k);
  xyz = M.xyz(no(:),1:3);
  [n,n3]=size(xyz);
  if n<2, 
    disp('Warning a spline with less then 2 points - not rendered');
  else
    xyz = cubicsplinecurve(xyz,n*5,'equal');
    
    if h & ishandle(h),
      [msize,colo,mark]=objproperties('curve',typ,0);
      set(h,'xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3),'linewidth',msize,'color',colo);
      set(h,'userdata',k,'tag','curve');
    else
      [msize,colo,mark]=objproperties('curve',typ,0);
      M.curvehandles(k) =...
	  line('xdata',xyz(:,1),'ydata',xyz(:,2),'zdata',xyz(:,3),'linewidth',msize,'color',colo);
      set(M.curvehandles(k),'userdata',k,'tag','curve');
    end;
    M=update_curvepoint(Ko,M,k);  % update points associate with the curve.
  end;
end;
% now check if a curve point is defined relative to this curve:
[ma,na]=size(M.A);
for k=1:ma,
  if M.edgetable(k,1) == 0,
    nv = compact(M.A,k);
    nv1 = abs(min(nv));
    nv2 = abs(max(nv));
    nn=M.vertextable([nv1,nv2]);
    X = [M.xyz(nn(1),1:3);M.xyz(nn(2),1:3)];
    h = M.edgehandles(k);
    if h~=0 & ishandle(h),
      set(h,'xdata',X(:,1),'ydata',X(:,2),'zdata',X(:,3));
    else
      M.edgehandles(k) = line('xdata',X(:,1),'ydata',X(:,2),'zdata',X(:,3),'color','g');
      set(M.edgehandles(k),'tag','straightedge');
      set(M.edgehandles(k),'hittest','off');         % initial 
    end;
  end;
end;
M.redpoint=0;
M.lastpoint=0;
M.redcurve=0;
return;  % rerender_graph
  
%-----------------------------------------------------------
function M= save_binary(Ko,M)
 Mobj = M;
 l = length(M.pointhandles); Mobj.pointhandles=zeros(l,1);
 l = length(M.curvehandles); Mobj.curvehandles=zeros(l,1);
 l = length(M.edgehandles); Mobj.edgehandles = zeros(l,1);
 save 'KObject' Mobj
return;
%-----------------------------------------------------------
function M=load_binary(Ko);
load 'KObject'
M = Mobj;
l = length(M.pointhandles); M.pointhandles=zeros(l,1);
l = length(M.curvehandles); M.curvehandles=zeros(l,1);
l = length(M.edgehandles); M.edgehandles = zeros(l,1);
M=rerender_graph(Ko,M);
M.freecell=0;         % gnadenlos....
M.freecellnumber=0;   % impitoyablement ...
M.macrocoords=[];     % merciless...  resetreset.
M.macrosplines = [0,0]; % 
M.specialcell=0;  % special cell: 1 cell treated as multi cell macro
M.autoupdate=1;
M.mouseoperation=0;
M.gluons=0;
M.lastcurve=0;   % consistent with later changes
M.lastpoint=0;   % consistent with later changes
M=clear_object_memory(Ko,M,'points');   
M=clear_object_memory(Ko,M,'curves');   
if ~isfield(M,'celllabels'),
  nc = length(M.cells);
  M.celllabels={};
  for k=1:nc,
    M.celllabels{k} = 'xxxxxxxx';
  end;
  tissueMaster('index','xxxxxxxx');
end;

return;
%-----------------------------------------------------------
function kill_surfaces(M);
a = findobj('tag','loopsurface');
for k=1:length(a),
  delete(a(k));
end;
return;
%-----------------------------------------------------------
function M=make_skinsurfaces(Ko,M)
sectionskin(findobj('tag','FEMKonstructorAxes'),0);  % external subroutine
return;
%--------------------------------------------------------------

function minimize_visible(M,how)
% if a bunch of cells get selected, checks which vertices and edges should be visible.
% how can be on or off.
  switch how,
   case 'on',
    he = findobj('tag','straightedge');
    if ~isempty(he), set(he,'visible','on'); end;
    he = findobj('tag','curve');
    if ~isempty(he), set(he,'visible','on'); end;
    he = M.pointhandles;
    if ~isempty(he), set(he,'visible','on'); end;
    
   case 'off',
    npts=length(M.pointhandles);
    markpts = zeros(npts,1);
    [nedge,n3]=size(M.edgetable);
    markedge = zeros(nedge,1);
    [markloops,surfcolor]=set_surface_colors(M);  % see below
    markloops = find(markloops==1); 
    
    for k=1:length(markloops),
      nsurf = markloops(k);
      [seqpts,seqedge] = LoopTopology(M,nsurf);     % sequence of verticies and edges
                                                    % with sign.
      seqpts=abs(seqpts);
      seqedge=abs(seqedge);
      markpts(seqpts)=1;
      markedge(seqedge) = 1;
    end;
    markpts=find(markpts==0);
    set(M.pointhandles(markpts),'visible','off');
    he = findobj('tag','curve');
    if ~isempty(he), set(he,'visible','off'); end;
    he = findobj('tag','straightedge');
    if ~isempty(he), set(he,'visible','off'); end;
    
    onedges = find(markedge==1);  % the ones that should stay on.
    hands = M.edgehandles(onedges);  % if they are straight lines
    set(hands(find(hands)),'visible','on');
    oncurvs = M.edgetable(find((markedge==1) & (M.edgetable(:,1) ~= 0)),1);   % the curves
    set(M.curvehandles(oncurvs),'visible','on');
    % finally switch on additional control points in the splines.
    for k=1:length(oncurvs),
      no = M.curves{oncurvs(k)};
      typ=M.pointtypes(no);
      nn = find(typ==3);
      if ~isempty(nn),
	no = no(nn);
	set(M.pointhandles(no),'visible','on');
      end;
    end;
  end;
  return;
    

function [mark,surfcolor]=set_surface_colors(M)
% in renderstruct one has a bunch of names and color names
% We search the cell labels for the names and attach a color to 
% each surface that is going to be drawn. 
  
  [mc,nc]=size(M.C);
  mark=zeros(nc,1); 
  surfcolor=cell(nc,1);
  
  colortab = colortableau;              % external subroutine that returns a
                                        % structure. 
  
  if isempty(colortab),
    for k=1:nc, 
      mark(k) = length(find(M.C(:,k)))==1; 
      surfcolor{k}='b';
    end;
  else,
    for k=1:length(M.cells),
      cle = M.celllabels{k};
      drawit=0;
      for j=1:length(colortab),
	if isequal(colortab{j}{1}, cle),           % compare names
	  col=colortab{j}{2};
	  drawit=tissueMaster('draw?',cle,k);       % returns 0 or 1.
	end;
      end;
      	  
      if drawit==1,
	llops = abs(compact(M.C,k));
	mark(llops) = mark(llops)+1;              % if marked once only, it will be drawn
	for j=1:length(llops),
	  surfcolor{llops(j)} = col;
	end;
      end;
    end;
  end;
  return;
  

function M=make_surfaces(Ko,M)
% for all loops that are accessed only once.
kill_surfaces(M);
%[mc,nc]=size(M.C);
%mark=zeros(nc,1);
%for k=1:nc, mark(k) = length(find(M.C(:,k)))==1; end;
% mark contains a 1 if the surface was only once mentioned in C, hence is
% an outer surface.
[mark,surfcolor]=set_surface_colors(M);

mark=find(mark==1);
% for all loops that are exterior (used only once), we calculate the 
% point and edge sequences. Then we construct a data structure that can
% be used by the coonspatch algorithm.
for k=1:length(mark),
  nsurf = mark(k);
  scolor = surfcolor{nsurf};
  
  [seqpts,seqedge] = LoopTopology(M,nsurf);     % sequence of verticies and edges with sign.
  nlin=length(seqedge);
  vonlin=zeros(nlin,2);
  for j=1:nlin,
    ed = seqedge(j);
    med=abs(ed);        % index into M.edgetable
    if ed>0,
      vonbis(j,1) = M.edgetable(med,2);
      vonbis(j,2) = M.edgetable(med,3);
    else
      vonbis(j,2) = M.edgetable(med,2);
      vonbis(j,1) = M.edgetable(med,3);
    end;    
    cuindx = M.edgetable(med,1);
    if cuindx==0,
      cuvs{j}=[];
    else
      cuvs{j} = M.curves{cuindx};
    end;
  end;
    
  nstuetz=set_stuetz(M);
  
  if nlin==3,
    [rx,ry,rz]=triangleface(M.xyz,cuvs,vonbis,nstuetz);
    ix = triangleindex(nstuetz);
    li=patch(rx(ix'),ry(ix'),rz(ix'),'b');
    set(li,'facecolor',scolor,'edgecolor','none');
    set(li,'tag','loopsurface','hittest','off');
  else,
    [rx,ry,rz]=bilinearcoonsplus(M.xyz,cuvs,vonbis,nstuetz);      % automatically nstuetzxnstuetz plates.
    li = surface(rx,ry,rz);
    set(li,'facecolor',scolor,'edgecolor','none','hittest','off');
    
%    li = line('xdata',rx(:),'ydata',ry(:),'zdata',rz(:),...
%    'marker','.','linestyle','none','markersize',1,'color','g');
    set(li,'tag','loopsurface','hittest','off');
  end;  
end;
%-----------------------------------------------------------
function M=light_operation(M,what,parms)
  if nargin<3, parms=[]; end;
  switch what,
   case 'on',
    o = findobj('tag','Werkstattlicht');
    if isempty(o),
      o(1) = light;
      o(2) = light('position',-get(o(1),'position'),'color',[0.7 0.5 0.5]);
      set(o,'tag','Werkstattlicht');
    end;
   case 'off',
    o = findobj('tag','Werkstattlicht');
    if ~isempty(o),
      for j=1:length(o),
	if  ishandle(o(j)),
	  delete(o(j));
	end;
      end;
    end;
   case 'move',
    o = findobj('tag','Werkstattlicht');
    if ~isempty(o),
      if length(parms)==3,
	set(o,'position',parms);
      else
	set(o,'position',-get(o,'position'));
      end;
    end;
  end;
  return;
  

%-----------------------------------------------------------
function M = convert_to_bezier(Ko,M)
% Takes the current graph and creates a data structure
% This subroutine is fairly long: It does essentially the following: 
% First, for each loop it generates the interior points. Their number depends on a fixed
% parameter, nipol, indicating how many points are on an edge. The interior points of each 
% surface are registered as type 4 points with a reference to the index of the surface,
% hence making it possible to find to which surface a point belongs. For each surface
% a little record that contains information of size, type (triangular or quadratic) and
% the indicies of the interior points is generated and stored in the array M.bezsurfs.
% Similarly for all edges, nipol interior points are calculated (by curvepiece) and
% registered as type 3 points with a reference to the edge number. For each edge a small
% array is stored that contains the indices of its control vertices. (point indicies are
% indicies into M.xyz and M.pointreferences. 
% 

[ma,na]=size(M.A);
if ma==0, 
  return;
end

[mb,nb]=size(M.B);
nipol=4;       % for now a constant: number of point on each edge.

% und jetzt: For each edge in A. fill in a new edge description bezedge

M.bezedges = cell(ma,1);   % empty cell
M.bezsurfs = cell(mb,1);

for nsurf=1:mb,       %    for all surfaces (internals and external)
  [seqpts,seqedge] = LoopTopology(M,nsurf);   % external subroutine
  nlin=length(seqedge);   % number of edges
  vonlin=zeros(nlin,2);    
  for j=1:nlin,
    ed = seqedge(j);
    med=abs(ed);          % index into M.edgetable
    if ed>0,
      vonbis(j,1) = M.edgetable(med,2);
      vonbis(j,2) = M.edgetable(med,3);
    else
      vonbis(j,2) = M.edgetable(med,2);
      vonbis(j,1) = M.edgetable(med,3);
    end;    
    cuindx = M.edgetable(med,1);
    if cuindx==0,
      cuvs{j}=[];
    else
      cuvs{j} = M.curves{cuindx};
    end;
  end;
  newcoords=[]; num_newcoords=0;
  if nlin==3,
    [rx,ry,rz]=triangleface(M.xyz,cuvs,vonbis,nipol);
    ix = triangleindex(nipol);
    [rand,inner]=triangle_rand(nipol);      % addresses of inner points and boarder  points (rand)
    newcoords=[rx(inner),ry(inner),rz(inner)];   % new inner points (coordinates)
    nin=length(inner);
    bezinfo={'t',nipol};
  elseif nlin==4,
    [rx,ry,rz]=bilinearcoonsplus(M.xyz,cuvs,vonbis,nipol);   % automatically 5x5 plates.
    rxi = rx(2:nipol-1,2:nipol-1);
    ryi = ry(2:nipol-1,2:nipol-1);
    rzi = rz(2:nipol-1,2:nipol-1);
    [rand,inner]=viereck_rand(nipol);                % attention! inner is 2-dimensional
    newcoords = [rxi(:),ryi(:),rzi(:)];
    nin=(nipol-2)^2;                                 % number of internal points.
    bezinfo = {'q',nipol};
  end;  
  nindicies=zeros(nin,1);
  for j=1:nin,
    [M,pointindx]= addobject(Ko,M,'point',newcoords(j,1:3),4,nsurf);  % type=4, reference = nsurf
    nindicies(j)=pointindx;
    set(M.pointhandles(pointindx),'hittest','off','visible','off');   % initially not selectable
  end;
  M.bezsurfs{nsurf}={bezinfo,nindicies};
end;
% %%  M.surfacexyz = newcoords;
% 
% And now for the curve points
newcoords =[]; num_newcoords=0;

for nedge=1:ma,            % for each edge
  ll = compact(M.A,nedge);
  nv1=abs(min(ll));           % vertex number at beginning of edge
  nv2=max(ll);                % and at end ( must be positive.)
  np1=M.vertextable(nv1);
  np2=M.vertextable(nv2);     % get the actual point indices
  
  cuindx = M.edgetable(nedge,1);  % The spline curve on which these vertices are - or null
  if cuindx==0,
    cu=[];
  else
    cu=M.curves{cuindx};
  end;
  xc = curvepiece(M.xyz,cu,np1,np2,nipol); % returns nipol points on the curve or on straight line
  newcoords = xc(2:nipol-1,:);  % internal points only
  nin=nipol-2;                 %    number of points on edge not counting corner verticies.
  % create control vertices with information of edge number and index within the edge
  nindicies=zeros(nin,1);
  for j=1:nin,
    [M,pointindx]= addobject(Ko,M,'point',newcoords(j,:),3,nedge);  % type=3, reference = nedge
    nindicies(j)=pointindx;
    set(M.pointhandles(pointindx),'hittest','off','visible','off');   % initially not selectable
  end;
  M.bezedges{nedge} = nindicies;
end;

return;


function M=render_all_bezier(Ko,M)
% brings up all bezier surfaces.
[mb,nb]=size(M.B);
[mc,nc]=size(M.C);
mark=zeros(nc,1);
for k=1:nc, mark(k) = length(find(M.C(:,k)))==1; end;

for nsurf=1:mb,   % number of loops
  pa = render_bezier(Ko,M,nsurf,mark(nsurf));
  set(pa,'userdata',nsurf,'buttondownfcn','FEMKonstructor(''hit'',''surface'',get(gco,''userdata''));');
  M.beziersurfacehandles(nsurf) = pa;
end;
M.bezierdefined=1;
M.lasthitbezier=0;
return;

%-----------------------------------------------------------
function M=rerender_all_bezier(Ko,M)
% brings up all bezier surfaces.
[mb,nb]=size(M.B);
[mc,nc]=size(M.C);
mark=zeros(nc,1);
for k=1:nc, mark(k) = length(find(M.C(:,k)))==1; end;

for nsurf=1:mb,   % number of loops
  pa = M.beziersurfacehandes(nsurf);
  if ishandle(pa),
     rerender_bezier(Ko,M,nsurf,mark(nsurf));
  end;
end;
return;
%--------------------------------------------------------

function pa=render_bezier(Ko,M,nsurf,outer)
% creates the Bezier surfaces and stores their handles.
if nargin < 4, outer=1; end;
XB = get_surface_controlpoints(M,nsurf);
su = M.bezsurfs{nsurf};
typ=su{1}{1};    % letter 't' or 'q'
nipol=su{1}{2};  % number of control vertices on edge
nstuetz=set_stuetz(M);
if outer==0, color='r'; elseif outer==1, color='b'; end;
switch typ,
case 'q',           % rectangular
  X=allbezier('array2d',XB(:,:,1),XB(:,:,2),XB(:,:,3),[0:(nstuetz-1)]/(nstuetz-1),[0:(nstuetz-1)]/(nstuetz-1));  
  pa = surface(X(:,:,1),X(:,:,2),X(:,:,3));
  set(pa,'facecolor',color,'edgecolor','none','tag','beziersurface');
case 't',           % triangular
  X=allbezier('triinterp',XB,nipol,nstuetz);
  [rst,ix] = triangledomain(nstuetz);  
  pa = patch(X{1},X{2},X{3},'b');
  set(pa,'facecolor',color,'edgecolor','none','tag','beziersurface');
end;
return;

%-------------------------------------------------------
function rerender_bezier(Ko,M,nsurf,marked)
if M.bezierdefined==0, return; end;
XB = get_surface_controlpoints(M,nsurf);
su = M.bezsurfs{nsurf};
typ=su{1}{1};    % letter 't' or 'q'
nipol=su{1}{2};  % number of control vertices on edge
nstuetz=set_stuetz(M);
pa = M.beziersurfacehandles(nsurf);
switch typ,
case 'q',           % rectangular
  X=allbezier('array2d',XB(:,:,1),XB(:,:,2),XB(:,:,3),[0:(nstuetz-1)]/(nstuetz-1),[0:(nstuetz-1)]/(nstuetz-1));  
  set(pa,'xdata',X(:,:,1),'ydata',X(:,:,2),'zdata',X(:,:,3));
case 't',           % triangular
  X=allbezier('triinterp',XB,nipol,nstuetz);
  [rst,ix] = triangledomain(nstuetz);  
  set(pa,'xdata',X{1},'ydata',X{2},'zdata',X{3});
end;

return;


%--------------------------------------------------------

function XE=get_edge_controlpoints(M,ndge)
% edge can be positive or negative, if negative invert the sequence.
nedge=abs(ndge);
ll = compact(M.A,nedge);
vtx1=abs(min(ll));  % beginning vertex
vtx2= max(ll);      % end vertex
x1 = M.xyz(M.vertextable(vtx1),1:3);  % coordinates of beginning point
x2 = M.xyz(M.vertextable(vtx2),1:3);  % coordinates of end point
indxseq = M.bezedges{nedge};
xedg = M.xyz(indxseq(:),1:3);
XE = [x1;xedg;x2];
if ndge<0, XE = flipud(XE); end;
return; 
%--------------------------------------------------------

function hE=get_edge_controlhandles(M,ndge)
% edge can be positive or negative, if negative invert the sequence.
nedge=abs(ndge);
ll = compact(M.A,nedge);
vtx1=abs(min(ll));  % beginning vertex
vtx2= max(ll);      % end vertex
h1 = M.pointhandles(M.vertextable(vtx1));  % coordinates of beginning point
h2 = M.pointhandles(M.vertextable(vtx2));  % coordinates of end point
indxseq = M.bezedges{nedge};
hedg = M.pointhandles(indxseq(:));
hE = [h1;hedg,h2];
if ndge<0, hE = flipud(hE); end;
return;  

%--------------------------------------------------------

function XB=get_surface_controlpoints(M,nsurf)
% retrieve the surface control points of a Bezier surface. tricky tricky mess
su = M.bezsurfs{nsurf}; 
typ=su{1}{1};            % letter 't' or 'q'
nipol = su{1}{2};
pindicies = su{2};       % interior point indicies into M.surfacexyz
% to retrieve the edges of the surface:
[seqpts,seqedge] = LoopTopology(M,nsurf);   % seqedge: edge sequence
nlin=length(seqedge);                       % number of edges
vonlin=zeros(nlin,2); 

switch typ,
 case 't',
 XB = zeros((nipol*(nipol+1))/2,3);
 [rand,inner] = triangle_rand(nipol);
 case 'q',
 XB = zeros(nipol*nipol,3);
 [rand,inner] = viereck_rand(nipol);
end;

% First, get the edge points.
for j=1:nlin,
  nl = seqedge(j);                         % can be positive or negative indicating direction
  nla=abs(nl);
  if nl<0,
    indxseq = flipud(M.bezedges{nla});    % interior point indices.
  else
    indxseq = M.bezedges{nla};
  end;  
%  disp(indxseq);
% The first vertex of each edge in the loop is in seqpts, here are it's coordinates:
  X1 = M.xyz(M.vertextable(seqpts(j)),1:3);
% So the first point with interior points (not including last point) is the following path:
  X = [X1; M.xyz(indxseq,1:3)];     
  bord = rand{j};
  XB(bord,1:3) = X;                         % fill in according to storage sequence
end;  
% finally - (good grief!) - for the inner points.
X = M.xyz(pindicies,1:3);         % their coordinates
XB(inner(:),1:3) = X;
if typ=='q',  XB = reshape(XB,[nipol,nipol,3]); end;
return;



%--------------------------------------------------------
function [rand,inner]=triangle_rand(n)
% address fields for triangles - used in FEMKonstructor.
switch n,
 case 3,
%    6
%   4 5
%  1 2 3
rand1 = [1 2 3]'; rand2=[3 4 5]'; rand3=[6 4 1]'; inner=[];

 case 4,
%     10
%    8  9
%   5  6  7 
% 1  2  3  4
rand1=[1 2 3 4]'; rand2=[4 7 9 10]'; rand3=[10 8 5 1]'; inner=[6]';

case 5
%        15
%      13 14
%    10 11 12
%   6  7  8  9  
% 1  2  3  4  5
rand1=[1 2 3 4 5 ]'; rand2=[5 9 12 14 15]'; rand3 = [15 13 10 6 1]'; inner=[7 8 11]';

case 6,
%          21
%        19 20
%      16 17 18 
%    12 13 14 15
%   7  8  9 10 11
% 1  2  3  4  5  6
rand1=[1:6]'; rand2=[6 11 15 18 20 21]'; rand3 = [21 19 16 12 7 1]'; inner=[8 9 10 13 14 17]';

case 7,
%            28
%          26 27
%        23 24 25
%      19 20 21 22 
%    14 15 16 17 18
%   8  9 10 11 12 13
% 1  2  3  4  5  6  7 
rand1=[1:7]'; rand2=[7 13 18 22 25 27 28]'; rand3=[28 26 23 19 14 8 1]';
inner=[9 10 11 12 15 16 17 20 21 24]';
otherwise,
  disp('wrong call for triangle_rand');
end;
% Exclude last point of each path:
rand{1}=rand1(1:n-1); rand{2}=rand2(1:n-1); rand{3}=rand3(1:n-1);
return;
%--------------------------------------------------------

function [rand,inner]=viereck_rand(n)
% addresses for boarders and internal points in a square.
% example:
%     1     6    11    16    21
%     2     7    12    17    22
%     3     8    13    18    23
%     4     9    14    19    24
%     5    10    15    20    25        
% Hence, rand1=[1 2 3 4], rand2=[5 10 15 20], rand3=[25 24 23 22], rand4=[21 16 11 6];
% The rest is stored to inner as 3x3 array.
rand{1}=[1:n-1];
rand{2}=[1:n-1]*n;
rand{3}=[(n*n):-1:(n*(n-1)+2)];
rand{4}=([n:-1:2]-1)*n+1;
inner=reshape(1:n*n,n,n);
inner=inner(2:n-1,2:n-1)';
return;

% If one clicks on a surface, it becomes transparent and its control verticies change
% color. They become selectable and all other surfaces are not selectable. If then one of 
% the surfaces has been chosen all the others will not react to mouse clicks. So we need an
% additional command, a mouse button that brings this machine to its initial state
% again. At any rate, after the control vertex of a surface become selectable their
% coordinates can be changed in the changer. Once the changing is done, the surface is
% recomputed and displayed as transparent again. 
% 
% A little change here:
% Let's continue the planning: Initially all surfaces are clickable. There are external
% ones and internal ones - which are held in another color. If one clicks on a surface
% with the middle mouse, it becomes transparent. If one clicks on one with the left mouse
% button it also becomes transparent but in addition it reveils its control points. All
% other surfaces then become unclickable (le seul methode de les délivrer est de appuyer
% un buton) so that it is possible to select a control point. Once the control point is
% moved, immediately the other visible surfaces become clickable and the one for which a
% control point was selected becomes nontransparent. Thus the idea here is: click on a
% surface, select a vertex, change vertex and see immediately the effects on the surface.
% The representation of the surfaces while they are transparent should also be changed.
% If there are too many surfaces, the mesh of edges is too dense. Replace it by a mesh
% that only shows the control vertices..??
% The selection type of the figure is 'normal'  left mouse button
%                                     'extend'  middle
%                                     'alt'     right
%                                     'open'    doubleclick
%-------------------------------------------------

% the following function cannot be used since matlab has a goddam bug. It stacks the
% surface and patch objects simply in the sequence of their generation and the last object 
% created is the most probable to be picked if several overlap. In other words their
% z-buffering technique is crap. Furthermore, I notcied that if the OpenGL rendering is
% enabled, object picking isn't possible at all! So I have to make some goofy
% workaround. The one I came up with consists in the following: I just click on a surface
% until it gets finally selected and use a button event to really select it. In other
% words:
% I click on a surface, Matlab finds SOME surface, probably not the one that I clicked on.
% I make that surface unclickable and change its color. When finally Matlab finds the
% surface that I want, I click on an extra button and woops it gets selected. 
% The one that did the work, in case matlab ever works right, is in Kjunk.m


function  M = surface_indexing(Ko,M,command,nsurf)
% This subroutine gets called if a Bezier surface is clicked by the mouse.
if M.bezierdefined==0, return; end;
switch command,
case 'pick',
  if M.lasthitbezier > 0 ,
    set(M.beziersurfacehandles(M.lasthitbezier),'edgecolor','none'); 
  end;
  surfh = M.beziersurfacehandles(nsurf);
  set(surfh,'edgecolor','w','hittest','off');        % make not hittable and the grid visible
  M.lasthitbezier = nsurf;
  Kmessage(sprintf('clicked on surface no %i',nsurf));
  
case 'hide',
  set(M.beziersurfacehandles(M.lasthitbezier),'visible','off'); 

case 'show',
  set(M.beziersurfacehandles(M.lasthitbezier),'visible','on'); 

case 'transparent',
  set(M.beziersurfacehandles(M.lasthitbezier),'visible','on','facecolor','none');
  
case 'expose',
  nsurf = M.lasthitbezier;  
  if nsurf==0, return; end;
  % make all others not hittable
  for k=1:length(M.bezsurfs),
    set(M.beziersurfacehandles(k),'hittest','off');
  end;
  % make this one invisible
  set(M.beziersurfacehandles(nsurf),'visible','off'); 
  % and expose control points: 
  hx=get_all_surfhandles(M,nsurf,1);
  for k=1:length(hx),
    set(hx(k),'hittest','on','visible','on','buttondownfcn','FEMKonstructor(''hit'',''point'',get(gco,''userdata''));'); 
  end;

case 'off',
  for k=1:length(M.bezsurfs),
    set(M.beziersurfacehandles(k),'hittest','off','visible','on','edgecolor','none');
  end;

case 'reset',
  nsurf=  M.lasthitbezier;
  M.lasthitbezier=0;
  % make all others not hittable
  if nsurf,
    hx=get_all_surfhandles(M,nsurf,1);  % not including edge verticies
    for k=1:length(hx),
      set(hx(k),'hittest','off','visible','off');
    end;
  end;
  disp('OK: surfaces to pickable again');
  for k=1:length(M.bezsurfs),
    set(M.beziersurfacehandles(k),'hittest','on','visible','on','edgecolor','none');
  end;
  
end;  
return;
%--------------------------------------------------------------
function M=bezier_select(Ko,M,code)  % from button callback

switch code,
case 'modify',
  M=surface_indexing(Ko,M,'expose');
case 'hide',
  M=surface_indexing(Ko,M,'hide');
case 'show',
  M=surface_indexing(Ko,M,'show');
case 'reset',
  M=surface_indexing(Ko,M,'reset');
end;

     

%-------------------------------------------------------

function hX=get_all_surfhandles(M,nsurf,sans)
% retrieve the surface control handles of a Bezier surface. tricky tricky mess
% and a simple rewriting of another subroutine
if nargin==2, sans=0; end;  % if equal to 1, avoid edgevertices.
su = M.bezsurfs{nsurf}; 
typ=su{1}{1};            % letter 't' or 'q'
nipol = su{1}{2};
pindicies = su{2};       % interior point indicies into M.surfacexyz
% to retrieve the edges of the surface:
[seqpts,seqedge] = LoopTopology(M,nsurf);   % seqedge: edge sequence
nlin=length(seqedge);                       % number of edges
vonlin=zeros(nlin,2); 

switch typ,
 case 't',
 hdl = zeros((nipol*(nipol+1))/2,1);
 [rand,inner] = triangle_rand(nipol);
 case 'q',
 hdl= zeros(nipol*nipol,1);
 [rand,inner] = viereck_rand(nipol);
end;

% First, get the edge points.
hX=[];
for j=1:nlin,
  nl = seqedge(j);                         % can be positive or negative indicating direction
  nla=abs(nl);
  if nl<0,
    indxseq = flipud(M.bezedges{nla});    % interior point indices.
  else
    indxseq = M.bezedges{nla};
  end;  
% The first vertex of each edge in the loop is in seqpts, here are it's coordinates:
  if sans==0,  hX = [hX; M.pointhandles(M.vertextable(seqpts(j)))]; end;
% So the first point with interior points (not including last point) is the following path:
  hX = [hX; M.pointhandles(indxseq(:))];     
end;  
% finally - (good grief!) - for the inner points.
hX = [hX; M.pointhandles(pindicies(:))];         % their coordinates
return;

function M=make_bezedges(Ko,M)
nstuetz=set_stuetz(M);
s = [0:nstuetz-1]/(nstuetz-1);
for k=1:length(M.bezedges),
  xe = get_edge_controlpoints(M,k);
  x = allbezier('interp',xe(:,1),xe(:,2),xe(:,3),s);
  M.bezieredgehandles(k) = line('xdata',x(:,1),'ydata',x(:,2),'zdata',x(:,3),'color','w');
  set(M.bezieredgehandles(k),'tag','bezieredge','hittest','off');
end;

%-------------------------------------------
function M=rerender_bezedge(Ko,M,nedge)
if M.bezierdefined==0, return; end;
nstuetz=set_stuetz(M);
s = [0:nstuetz-1]/(nstuetz-1);
xe = get_edge_controlpoints(M,nedge);
x = allbezier('interp',xe(:,1),xe(:,2),xe(:,3),s);
set(M.bezieredgehandles(nedge),'xdata',x(:,1),'ydata',x(:,2),'zdata',x(:,3));
nsurfs = find(M.B(:,nedge));
for k=1:length(nsurfs),
  rerender_bezier(Ko,M,nsurfs(k));
end;
return;

function M=rerender_bezedges(Ko,M)
% actually finds the surfaces associate with the edge and rerenders those
if M.bezierdefined==0, return; end;
for k=1:length(M.bezedges),
  M=rerender_bezedge(Ko,M,k);
end;
return;


%-------------------------------------------
function M=rerender_bezEcke(Ko,M,necke)  
% find all surfaces associated with a vertex and rerender them
% 
if M.bezierdefined==0, return; end;
nlins = find(M.A(:,necke));
for k=1:length(nlins),
  M=rerender_bezedge(Ko,M,nlins(k));
end;
% next rerender all surfaces associated with the Eckenvertex.
nsurfs = surfaces_from_ecke(M,necke);
for k=1:length(nsurfs),
  rerender_bezier(Ko,M,nsurfs(k));
end;

return;

%---------------------------------------------
function nsurfs=surfaces_from_ecke(M,necke)
nlins=find(M.A(:,necke));
nsurfs=[];
for k=1:length(nlins),
  nsu=find(M.B(:,nlins(k)));
  nsurfs=[nsurfs;nsu];
end;
nsurfs = minimalset(nsurfs);
return;


%-------------------------------------------
function M=bezier_rendermode(Ko,M,nsurf,how)
% rendering a surface either transparent or opaque
h = M.beziersurfacehandles(nsurf);

switch how,
case 'opaqueON'
 set(h,'facecolor','edgecolor','none','b','hittest','on');
case 'opaquenOFF',
 set(h,'facecolor','edgecolor','none','b','hittest','off');
case 'transparent',
 set(h,'facecolor','none','edgecolor','w','hittest','off');
end;
return;

%-------------------------------------------
function toggle_oldlines(arg)
fo = findobj('tag','straightedge');
set(fo,'visible',arg);
fo = findobj('tag','curve');
set(fo,'visible',arg);
return;


%-------------------------------------------------
function M=translate_all(Ko,M,vec)
M.xyz(:,1) = M.xyz(:,1)+vec(1);
M.xyz(:,2) = M.xyz(:,2)+vec(2);
M.xyz(:,3) = M.xyz(:,3)+vec(3);
ax = findobj('tag','FEMKonstructorAxes');
l = findobj(get(ax,'children'),'type','line');
if ~isempty(l), delete_lines(l); end;
M = rerender_graph(Ko,M);
return;
%-------------------------------------------------
function M=rotate_all(Ko,M,R)
x = M.xyz(:,1:3);
[n,n3]=size(x);
mx=sum(x)/n;
mx=ones(n,1)*mx;
x = (x-mx)*R+mx;
M.xyz(:,1:3) = x;
ax = findobj('tag','FEMKonstructorAxes');
l = findobj(get(ax,'children'),'type','line');
if ~isempty(l), delete_lines(l); end;
M = rerender_graph(Ko,M);
return;
%-------------------------------------------------
function M=scale_all(Ko,M,v)
x = M.xyz(:,1:3);
[n,n3]=size(x);
mx=sum(x)/n;
x(:,1) = (x(:,1)-mx(1))*v(1) + mx(1);
x(:,2) = (x(:,2)-mx(2))*v(2) + mx(2);
x(:,3) = (x(:,3)-mx(3))*v(3) + mx(3);
M.xyz(:,1:3)=x;
ax = findobj('tag','FEMKonstructorAxes');
l = findobj(get(ax,'children'),'type','line');
if ~isempty(l), delete_lines(l); end;
M = rerender_graph(Ko,M);
return;
%-------------------------------------------------
function delete_lines(l)
% delete all lines that don't have the tag 'ebene'
for k=1:length(l),
  t = get(l(k),'tag');
  if ~isequal(t,'ebene'), 
    delete(l(k));
  end;
end;
return;

%------------------------------------------------
function c= make_a_label(name,n)
% assuming that there is only one '.'
if isempty(name), c = blanks(n); return; end;

np = find(name=='.');
if isempty(np), 
  c = name; 
else 
  c = name(1:np(1)-1) ;
end;
nl=length(c);
if nl>=n; 
  c=c(1:n); 
else
  c = [c,blanks(n-nl)];
end;
return;

%-------------------------------------------------

function [M,ret]=include_macro(Ko,M,filename,cellname)
fid = fopen(filename,'r');
if fid<0, Kmessage(['Desolé, le fichier ''''',filename,''''' n''existe pas!']); return; end;
ret=1;
eof=0;

if isempty(cellname),
  cellname=make_a_label(filename,8);   % create a  8 character  label
end;

lastspline=length(M.curvehandles);
nosplines=0;
nptprior = length(M.pointhandles);
vertprior=length(M.vertextable);
numvert=0;  % new vertex counter
previouslastcell=length(M.cells);
M.freecellnumber=0;  % counted up in add_a_cell.
coords=[];

normload=1;
fixedcell=0;   % 

while ~eof,
  [Car,n] = fscanf(fid,'%c',1);
  switch Car,
  case 'E',
   eof=1;
   break;
   
   case 'Y',            % cludge for special files: first line starts
                        % with letter Y
    normload=0;         %  
   case 'F',           
    M.specialcell=1;    % oh oh well. what a cludge 1 cell treated as macro.
    
   case 'V',            % a vertex with coordinates xyz
    [no,n]=fscanf(fid,'%i',1);
    [xyz,n]=fscanf(fid,'%f',3);
    [M,np]= addobject(Ko,M,'point',xyz(:)',1);   % is going to be a vertex.
    coords=[coords;xyz(:)'];                     % collect the coordinates to store extra
    numvert=numvert+1;
    M.vertextable = [M.vertextable,np];
    macrovertx(numvert)=np;   % update transform table. 
    
   case 'I',            % a for interpolation
    [no,n]=fscanf(fid,'%i',1);
    [xyz,n]=fscanf(fid,'%f',3);
    [M,np]= addobject(Ko,M,'point',xyz(:)',3);   % is going to be a vertex.
    coords=[coords;xyz(:)'];                     % collect the coordinates to store extra
    
   case 'S',            % a spline curve 
    [no,n] = fscanf(fid,'%i',1);   % find number of points.
    [no,n] = fscanf(fid,'%i',no);
			 %  disp(no');
			 no = no+nptprior;         % offset point numbers by already existing point indices.
			 [M,ret]=addobject(Ko,M,'curve',no,1); 
			 nosplines=nosplines+1;
			 
   case 'C',            % a cube.
    [no,n] = fscanf(fid,'%i',8);     
    if normload,
      no = no+vertprior;
    else
      for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
    end;
    M = add_a_cell(Ko,M,'cube',no',cellname);
    
   case 'P',            % a prism
    [no,n] = fscanf(fid,'%i',6);      
    if normload,
      no = no+vertprior;
    else
      for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
    end;
    M = add_a_cell(Ko,M,'prism',no',cellname);
    
   case 'T',            % a tetrahedron.
    [no,n] = fscanf(fid,'%i',4); 
    if normload,
      no = no+vertprior;
    else
      for ii=1:length(no), no(ii) = find(M.vertextable==no(ii)); end;
    end;
    M = add_a_cell(Ko,M,'tetra',no',cellname);
    
   case 'L',                        % cell label
    [no,n] = fscanf(fid,' %4i',1);
    [name,n]=fscanf(fid,'%8s',1);
    name = make_a_label(name,8);
    M.celllabels{previouslastcell+no} = name;
    tissueMaster('index',name);
  end;
end;
ncells=length(M.cells);
M.freecell = previouslastcell+1;   %  now we are dealing with multiple free cells.
fclose(fid);
setedit_values(Ko.topoedits,[0 0 0 0 0 0 10 10 10]);  % so the thing can be scaled and
                                                      % turned
						      % see move_free_cell
						      
M.macrocoords=[];
M.macrosplines = [lastspline,nosplines];
if ncells==0, return; end;
if normload == 1, [M,flg]=move_free_cell(Ko,M); end;  % otherwise it will
                                                      % translate it to the
                                                      % origin. In
                                                      % general, just one cell
                                                      % is considered to
                                                      % be a translatable cell.
						      % 
return;

%-------------------------------------------------
function M=add_a_cell(Ko,M,type,vindex,name)
[nb,mb]=size(M.B);
[M.A,M.B,M.C,newedges]=Glue(M.A,M.B,M.C, type, 0, vindex);  % extern (makes primitives).
[nbn,mbn]=size(M.B);
newloops = [(nb+1):1:nbn];    
% (loops are added at the end of the matrix).
% So we just got the new topology. Now we still need to create the edges.
% newedges  contains the row indices that where added to the A matrix.
% So we create a different representation of that:

% a cell array
% this contains the edges that need to be drawn. Since they have not been
% in the topology before (Glue returns only new lines), they all need to 
% be plotted, from the initial point to the end. 

if ~isempty(newedges), 
  tab = compact(M.A,newedges);  
  [M,flg]=make_edges(Ko,M,tab);
end;
% this first of all just draws curves, but gives them a special index. 
if ~isempty(newloops),
  [M,flg]=make_loops(Ko,M,newloops);
end;
ncell=length(M.cells);
M.cells{ncell+1} = {type, vindex};
M.celllabels{ncell+1} = name;
M.freecellnumber = M.freecellnumber+1;
return;


function M=make_a_cludge(Ko,M)
% This function should usually not be used but I bet it's going to be used all the time to 
% make some ugly changes. The first time it was used was to convert all control vertices
% that are not in vertextable into points of type 3.
for k=1:length(M.pointhandles),
  ni = find(M.vertextable==k);
  if isempty(ni),
    M.pointtypes(k)=3;
  end;
end;
return;

function Ko= select_view(Ko,M,code)
% affects only the viewing buttons.

if code==0,  % store the current view to the last one of the three viewing buttons if 
             % that had been pressed before. 
  v = get(gca,'view');
  nv=Ko.currentview;
  if nv>0,
    set(Ko.viewbuttons(nv),'userdata',v);
    set(Ko.viewbuttons(nv),'backgroundcolor','m');
  end;
elseif code == 1 | code == 2 | code ==3,
  v = get(Ko.viewbuttons(code),'userdata');
  if ~isempty(v),
    set(gca,'view',v);
    for j=1:3,
      u = get(Ko.viewbuttons(j),'userdata');
      if ~isempty(u),
	if j==code, co='m'; else, co='y'; end;
      else
	co='w';
      end;
      set(Ko.viewbuttons(j),'backgroundcolor',co);
    end;
  end;
  Ko.currentview=code;
end;


function M=activate_near_points(Ko,M,maxdist)
  dis = ebenenmaster('distance',M.xyz);
  dis = abs(dis)-maxdist;
  n = find(dis<0);
  if ~isempty(n),
    set(M.pointhandles,'hittest','off','visible','off');
    for k=1:length(n),
      set(M.pointhandles(n(k)),'hittest','on','color','r','visible','on');
    end;
  end;
  return;

function n = set_stuetz(M)
 n = 40;
 return;