blob = h;

blob_center = (size(blob) + 1) / 2


T1 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    -blob_center 1];

theta = pi/8;

T2 = [cos(theta)  0      -sin(theta)   0
    0             1              0     0
    sin(theta)    0       cos(theta)   0
    0             0              0     1]

T3 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    blob_center 1]

T = T1 * T2 * T3

tform = maketform('affine', T);

tformfwd(blob_center, tform)

R = makeresampler('linear', 'fill');

TDIMS_A = [1 2 3];
TDIMS_B = [1 2 3];
TSIZE_B = size(blob);
TMAP_B = [];
F = 0;
blob2 = tformarray(blob, tform, R, TDIMS_A, TDIMS_B, TSIZE_B, TMAP_B, F);
hs2 = tformarray(hs, tform, R, TDIMS_A, TDIMS_B, TSIZE_B, TMAP_B, F);

h = blob2;

% First set all function values not close to surface far away from 0
disp('find'); tic;
h(find(abs(h)>.01))  = 2;
toc;

% Find all points on the surface lying outside a ball of radius r
disp('intersect'); tic;
h(intersect(find(h < 2),find(hs2 >= .51))) = 2;
toc;


clf
isosurface(h,.1)
%p = patch(isosurface(blob2,0));
%set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
%view(3)
camlight
%lighting gouraud
