Index: Damkjer/ImageProcessing/BB_Inpainting/bb_diffuse.m
===================================================================
--- Damkjer/ImageProcessing/BB_Inpainting/bb_diffuse.m	(revision 0)
+++ Damkjer/ImageProcessing/BB_Inpainting/bb_diffuse.m	(revision 4)
@@ -1,23 +1,24 @@
-function O=bb_diffuse(I,R,k,t,idxs)
+function O=bb_diffuse(I,R,k,t,idxs,mask)
 
 H=fspecial('gaussian');
 I=imfilter(I,H,'replicate');
-[Ix, Iy] = grad(I);
-[Ixx,Ixy,Iyy] = fdiff2(I);
-Ix2 = Ix.^2;
-Iy2 = Iy.^2;
-norm2 = Ix2+Iy2+1e-10;
 
-g=(1./(1+k.*norm2));
-
-num=(Ix2.*Iyy-2.*Ix.*Iy.*Ixy+Iy2.*Ixx);
-sub=(1-g).*(I-R);
-
-diff = (g).*(num./norm2)-sub;
-
-O=min(max(0,I+t.*diff),1);
-
+   [Ix, Iy] = grad(I);
+   [Ixx,Ixy,Iyy] = fdiff2(I);
+   Ix2 = Ix.^2;
+   Iy2 = Iy.^2;
+   norm2 = Ix2+Iy2+1e-10;
+   
+   g=(1./(1+k.*norm2));
+   
+   num=(Ix2.*Iyy-2.*Ix.*Iy.*Ixy+Iy2.*Ixx);
+   sub=(1-g).*(I-R);
+   
+   diff = (g).*(num./norm2)-sub;
+   
+   O=min(max(0,I+t.*diff),1);
+   
 if (nargin<5)
-    return
+   return
 end
 
@@ -27,4 +28,10 @@
 Iy2 = Iy.^2;
 norm2 = Ix2 + Iy2 + 1e-10;
+
+O=I;
+
+%diff=Ix2.*Iyy-2.*Ix.*Iy.*Ixy+Iy2.*Ixx./norm2;
+%O(mask)=I(mask) + t * diff;
+
 for idx=1:size(idxs,1)
     x=idxs(idx,2);
Index: Damkjer/ImageProcessing/BB_Inpainting/bb_inpaint.m
===================================================================
--- Damkjer/ImageProcessing/BB_Inpainting/bb_inpaint.m	(revision 0)
+++ Damkjer/ImageProcessing/BB_Inpainting/bb_inpaint.m	(revision 4)
@@ -5,6 +5,6 @@
 % Parse inputs
 p = inputParser;
-defaultPaintSteps = 10;
-defaultSmoothSteps = 10;
+defaultPaintSteps = 15;%10
+defaultSmoothSteps = 2;%10
 defaultPaintTimeSlice = 0.02;
 defaultSmoothTimeSlice = 0.2;
@@ -32,29 +32,42 @@
 %Determine Omega (the inpaint region)
 mask=zeros(size(I));
-mask(:)=isnan(I(:));
+mask=isnan(I);
 
-P=I/255;
-P(isnan(I))=0.5;
+%P=I/255;
+%P=I;
+%P(mask)=0.5;
+P=bscb_preproc(I);
 
 J=P;
 R=P;
 
-[r,c]=ind2sub(size(I),find(mask));
-mask=[r,c];
+[r,c]=ind2sub(size(I),find(mask(:)));
+mask2=[r,c];
+
+size(mask2)
+
+h = waitbar(0,'Pre-Processing...');
 
 % Pre-processing: the whole image undergoes anisotropic diffusion smoothing
 P=bb_diffuse(P,R,k,ts);
 
+close(h);
+
 h = waitbar(0,'In Painting...');
 
 for pass=1:T;
-    for i=1:A
-        % Equation (8)
-        [Ix,Iy] = grad(P, mask);
-        [Ixxx,Ixxy,Ixyy,Iyyy]=fdiff3(P,mask);
+
+   for i=1:A
+      % Equation (8)
+        [Ix,Iy] = grad(P, mask2);
+        [Ixxx,Ixxy,Ixyy,Iyyy]=fdiff3(P,mask2);
+
+%        It=Ix.*(Ixxy+Iyyy)-Iy.*(Ixxx+Ixyy);
+%        P(mask)=min(max(0,P(mask)+ti.*It),1);
+
         
-        for idx=1:size(mask,1)
-            x=mask(idx,2);
-            y=mask(idx,1);
+        for idx=1:size(mask2,1)
+            x=mask2(idx,2);
+            y=mask2(idx,1);
 
             It=Ix(idx)*(Ixxy(idx)+Iyyy(idx))-Iy(idx)*(Ixxx(idx)+Ixyy(idx));
@@ -65,8 +78,8 @@
     
     for j = 1:B
-        P=bb_diffuse(P,R,k,ts,mask);
+       P=bb_diffuse(P,R,k,ts,mask2,mask);
     end
     
-    J(isnan(I))=P(isnan(I));
+    J(mask)=P(mask);
     
     if (mod(pass,10) == 0||pass==1||pass==T)
Index: Damkjer/ImageProcessing/BSCB_Inpainting/bscb_inpaint.m
===================================================================
--- Damkjer/ImageProcessing/BSCB_Inpainting/bscb_inpaint.m	(revision 0)
+++ Damkjer/ImageProcessing/BSCB_Inpainting/bscb_inpaint.m	(revision 4)
@@ -1,5 +1,16 @@
+% [OUT]=BSCB_INPAINT(IN,...) Bartalmio et al image inpainting.
+%   Performs image in-painting as described by Bertlamio, Sapiro, Caselles, and 
+%   Ballester (2000).
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2012-AUG-29   K. Damkjer
+%         Initial Coding.
+%      2013-FEB-04   K. Damkjer
+%         Updated to remove mis-interprestation of paper and add pre-processing
+%         step.
+%
 function J=bscb_inpaint(I,varargin)
-% BSCB_INPAINT: in-paints over nans in an matrix
-% usage: B=BSCB_INPAINT(A)          % default method
 
 % Parse inputs
@@ -7,6 +18,6 @@
 defaultPaintSteps = 15;
 defaultSmoothSteps = 2;
-defaultTimeSlice = 0.1;
-defaultIters = 3000;
+defaultTimeSlice = 0.2;
+defaultIters = 300;
 
 addRequired(p, 'I', @isfloat);
@@ -28,6 +39,8 @@
 mask(:)=isnan(I(:));
 
-P=I/255;
-P(isnan(I))=0.5;
+%P=I;
+%P(isnan(I))=0.5;
+
+P=bscb_preproc(I);
 
 J=P;
@@ -96,10 +109,10 @@
                 It = beta*sqrt(mxb^2+mxf^2+myb^2+myf^2);
 
-                mxb = max(0,Ixb(idx));
-                mxf = min(0,Ixf(idx));
-                myb = max(0,Iyb(idx));
-                myf = min(0,Iyf(idx));
+%                mxb = max(0,Ixb(idx));
+%                mxf = min(0,Ixf(idx));
+%                myb = max(0,Iyb(idx));
+%                myf = min(0,Iyf(idx));
                 
-                It = 0.5*(It + (-beta)*sqrt(mxb^2+mxf^2+myb^2+myf^2));
+%                It = 0.5*(It + (-beta)*sqrt(mxb^2+mxf^2+myb^2+myf^2));
             else
                 mxb = max(0,Ixb(idx));
@@ -110,10 +123,10 @@
                 It = beta*sqrt(mxb^2+mxf^2+myb^2+myf^2);
 
-                mxb = min(0,Ixb(idx));
-                mxf = max(0,Ixf(idx));
-                myb = min(0,Iyb(idx));
-                myf = max(0,Iyf(idx));
+%                mxb = min(0,Ixb(idx));
+%                mxf = max(0,Ixf(idx));
+%                myb = min(0,Iyb(idx));
+%                myf = max(0,Iyf(idx));
 
-                It = 0.5*(It + (-beta)*sqrt(mxb^2+mxf^2+myb^2+myf^2));
+%                It = 0.5*(It + (-beta)*sqrt(mxb^2+mxf^2+myb^2+myf^2));
 
             end
@@ -133,5 +146,4 @@
 %             end
             
-%             P(y,x)=min(max(0,P(y,x)+t*si*sqrt(si*It)),1);
             P(y,x)=min(max(0,P(y,x)+t*si*sqrt(sqrt(si*It))),1);
             
@@ -139,5 +151,5 @@
             % often it failed to propagate information into the inpainting
             % domain.
-%             P(y,x)=min(max(0,P(y,x)+t*0.5*It),1);
+%             P(y,x)=min(max(0,P(y,x)+t*It),1);
         end
     end
Index: Damkjer/ImageProcessing/BSCB_Inpainting/bscb_preproc.m
===================================================================
--- Damkjer/ImageProcessing/BSCB_Inpainting/bscb_preproc.m	(revision 4)
+++ Damkjer/ImageProcessing/BSCB_Inpainting/bscb_preproc.m	(revision 4)
@@ -0,0 +1,35 @@
+% [OUT]=BSCB_PREPROC(IN,...) Pre-processing for BSCB in-painting.
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2013-FEB-04   K. Damkjer
+%         Initial Coding.
+%
+function J=bscb_preproc(Image,varargin)
+
+% Parse inputs
+p = inputParser;
+
+addRequired(p, 'Image', @isfloat);
+
+% Establish default parameters
+parse(p,Image,varargin{:});
+Image=p.Results.Image;
+
+J=Image;
+
+%Determine Omega (the inpaint region)
+mask=zeros(size(J));
+mask(isnan(J))=1;
+labels=bwlabel(mask,8);
+
+regions=max(max(labels));
+
+for region=1:regions
+   roi=(labels==region);
+   del=bwmorph(roi,'dilate',1);
+   avg=nanmean(Image(del==1));
+   std=nanstd(Image(del==1));
+   J(roi==1)=avg+std.*randn(size(J(roi==1)));
+end
Index: Damkjer/ImageProcessing/Common/fdiff2.m
===================================================================
--- Damkjer/ImageProcessing/Common/fdiff2.m	(revision 0)
+++ Damkjer/ImageProcessing/Common/fdiff2.m	(revision 4)
@@ -1,5 +1,20 @@
+%FDIFF3 Numerical Second Differential.
+%   [FXX,FXY,FYY] = FDIFF2(F) returns the numerical second partial derivative of
+%   the matrix F. The spacing between points in each direction is assumed to be
+%   one.
+%
+%   [FXX,FXY,FYY] = FDIFF3(F,IDXS), where IDXS is a 2-D array, computes
+%   numerical third partial derivative at the locations specified by IDXS.
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2012-AUG-29   K. Damkjer
+%         Initial Coding.
+%      2013-FEB-04   K. Damkjer
+%         Updated output sizes or indexed case.
+%
+
 function [Oxx,Oxy,Oyy] = fdiff2( I,idxs )
-%FDIFF2 Summary of this function goes here
-%   Detailed explanation goes here
 
 if (nargin<2)
@@ -34,7 +49,7 @@
     Oyy=(Icp-2*I+Icn);
 else
-    Oxx = zeros(size(idxs));
-    Oxy = zeros(size(idxs));
-    Oyy = zeros(size(idxs));
+    Oxx = zeros(size(idxs,1),1);
+    Oxy = zeros(size(idxs,1),1);
+    Oyy = zeros(size(idxs,1),1);
 
     for idx=1:size(idxs,1)
Index: Damkjer/ImageProcessing/Common/fdiff3.m
===================================================================
--- Damkjer/ImageProcessing/Common/fdiff3.m	(revision 0)
+++ Damkjer/ImageProcessing/Common/fdiff3.m	(revision 4)
@@ -1,5 +1,20 @@
+%FDIFF3 Numerical Third Partial.
+%   [FXXX,FXXY,FXYY,FYYY] = FDIFF3(F) returns the numerical third partial
+%   derivative of the matrix F. The spacing between points in each direction is
+%   assumed to be one.
+%
+%   [FXXX,FXXY,FXYY,FYYY] = FDIFF3(F,IDXS), where IDXS is a 2-D array, computes
+%   numerical third partial derivative at the locations specified by IDXS.
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2012-AUG-29   K. Damkjer
+%         Initial Coding.
+%      2013-FEB-04   K. Damkjer
+%         Updated output sizes or indexed case.
+%
+
 function [Oxxx,Oxxy,Oxyy,Oyyy] = fdiff3( I,idxs )
-%FDIFF3 Summary of this function goes here
-%   Detailed explanation goes here
 
 if (nargin<2)
@@ -36,8 +51,8 @@
     Oyyy=0.5.*(-IcP+2.*Icp-2.*Icn+IcN);
 else
-    Oxxx = zeros(size(idxs));
-    Oxxy = zeros(size(idxs));
-    Oxyy = zeros(size(idxs));
-    Oyyy = zeros(size(idxs));
+    Oxxx = zeros(size(idxs,1),1);
+    Oxxy = zeros(size(idxs,1),1);
+    Oxyy = zeros(size(idxs,1),1);
+    Oyyy = zeros(size(idxs,1),1);
 
     for idx=1:size(idxs,1)
Index: Damkjer/ImageProcessing/Common/grad.m
===================================================================
--- Damkjer/ImageProcessing/Common/grad.m	(revision 0)
+++ Damkjer/ImageProcessing/Common/grad.m	(revision 4)
@@ -1,4 +1,3 @@
-function [Ox,Oy]=grad(I,idxs)
-%GRAD Gradient.
+%GRAD Numerical Gradient.
 %   [FX,FY] = GRAD(F) returns the numerical central gradient of the
 %   matrix F. FX corresponds to dF/dx, the differences in x (horizontal) 
@@ -13,12 +12,21 @@
 %   dimension of F, going across columns.  The second output FY is always
 %   the gradient along the 1st dimension of F, going across rows.
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2012-AUG-29   K. Damkjer
+%         Initial Coding.
+%      2013-FEB-04   K. Damkjer
+%         Updated output sizes or indexed case.
+%
 
-%   Copyright 2012 Kristian L. Damkjer
+function [Ox,Oy]=grad(I,idxs)
 
 if (nargin<2)
     [Ox,Oy]=gradient(I);
 else
-    Ox=zeros(size(idxs));
-    Oy=zeros(size(idxs));
+    Ox=zeros(size(idxs,1),1);
+    Oy=zeros(size(idxs,1),1);
 
     for idx=1:size(idxs,1)
Index: Damkjer/ImageProcessing/Common/isophote.m
===================================================================
--- Damkjer/ImageProcessing/Common/isophote.m	(revision 0)
+++ Damkjer/ImageProcessing/Common/isophote.m	(revision 4)
@@ -1,2 +1,17 @@
+%ISOPHOTE Compute isophote directions.
+%   [NX,NY] = ISOPHOTE(F) returns the estimate of the isophote direction as the
+%   direction perpendicular to gradient. Orientation arbitrarily selected as CW
+%   rotation of gradient (alternative would be CCW rotation of gradient).
+%
+%   Copyright 2012 Kristian L. Damkjer
+%
+%   Software History:
+%      2012-AUG-29   K. Damkjer
+%         Initial Coding.
+%      2013-FEB-04   K. Damkjer
+%         Updated to remove normalization to improve convergence and stability.
+%         Should add an option to enable and disable normalizing.
+%
+
 function [Ox,Oy]=isophote(Ix,Iy)
 
@@ -5,5 +20,6 @@
 
 for idx=1:length(Ix)
-    M=sqrt(Ix(idx)^2+Iy(idx)^2+1e-6);
+%    M=sqrt(Ix(idx)^2+Iy(idx)^2+1e-6);
+    M=1;
 
     Ox(idx) = -Iy(idx) / M;