source_dir = 'Sep 29 2003/';
names_like = 'dscn0*.jpg';
width = 2 + 7/32;
height = 2 + 3/32;
sep = 0.02;
xy_reduce = 5;

load('Sep 29 2003/outside.mat');
load('Sep 29 2003/med_tooth.mat');
load('Sep 29 2003/XII_VIII.mat');
load('Sep 29 2003/XIV.mat');

fig_3d = figure;
fig_pic = figure;
set(fig_pic, 'XDisplay', 'localhost:0.1');

files = dir([source_dir names_like]);
num_files = size(files, 1);

current_jpeg = imread([source_dir files(1).name]);
image(current_jpeg);
axis image;

title('Click at the top of the block.');
[dummy top] = ginput(1); 
title('Click at the bottom of the block.');
[dummy bottom] = ginput(1); 
title('Click at the left of the block.');
[left dummy] = ginput(1); 
title('Click at the right of the block.');
[right dummy] = ginput(1); 


V = zeros(round(size(current_jpeg, 2) / xy_reduce), round(size(current_jpeg, 1) / xy_reduce), num_files);
vert = zeros(size(V));
horz = zeros(size(V));

for file = 1 : num_files,
  current_jpeg = imread([source_dir files(file).name]);
  figure(fig_pic);
  clf;
  image(current_jpeg);
  title('Click left to add point, right to finish spline, middle for next picture.');
  axis image;
  hold on;
  adding_picture = true;
  while adding_picture,
    Xtemp = [];
    Ytemp = [];    
    adding_curve = true;
    while adding_curve,
      [y, x, button] = ginput(1);
      if button == 2, 
        adding_picture = false;
        adding_curve = false;
      elseif button == 3,
        adding_curve = false;
        Xtemp = [Xtemp(end), Xtemp, Xtemp(1)];
        Ytemp = [Ytemp(end), Ytemp, Ytemp(1)];
        t = 1 : size(Xtemp, 2);
        ts = 1 : 0.0005 : size(Xtemp, 2);
        xys = spline(t, [Xtemp; Ytemp], ts);
        
        x = round(xys(1,:) / xy_reduce);
        y = round(xys(2,:) / xy_reduce);
        for i = 1 : size(xys, 2),
          V(y(i), x(i), file) = 1;
        end;
        
        for i = 1 : size(V, 2),
          I = find(V(:, i, file));
          if any(I);
            vert(min(I) : max(I), i, file) = 1;
          end;
        end;
        
        for i = 1 : size(V, 1),
          I = find(V(i, :, file));
          if any(I);
            horz(i, min(I) : max(I), file) = 1;
          end;
        end;
        
        plot(xys(2,:), xys(1,:), 'b-');
      elseif button == 1,      
        Xtemp = [Xtemp, x];
        Ytemp = [Ytemp, y];
        plot(y, x, 'ro');
      end;
    end;
  end;
  figure(fig_3d);
  try,
    V = V | (vert & horz);
    clf;
    fv = isosurface(V, 0.5, 'verbose');
    fvi = fv;
    fvi.vertices(:, 1) = fvi.vertices(:, 1) * xy_reduce * width / (right - left);
    fvi.vertices(:, 2) = fvi.vertices(:, 2) * xy_reduce * height / (bottom - top);
    fvi.vertices(:, 3) = fvi.vertices(:, 3) * sep;
    patch(fvi_out, 'FaceLighting','Phong','EdgeColor','none','FaceColor','b','FaceAlpha',0.2);
    patch(fvi_med_tooth, 'FaceLighting','Phong','EdgeColor','none','FaceColor','g');
    patch(fvi_XII_VIII, 'FaceLighting','Phong','EdgeColor','none','FaceColor','c');
    patch(fvi_XIV, 'FaceLighting','Phong','EdgeColor','none','FaceColor','m');
    patch(fvi, 'FaceLighting','Phong','EdgeColor','none','FaceColor','r');
    axis equal
    light
    view([90 15]);
  catch,
    lasterr
  end;
  
end;