Spatial Bias on Interpolator (Doc & test & correction)

MMVII/Doc/Programmer/Interpolators.tex

@@ -122,7 +122,7 @@ \subsection{$1d$ interpolation kernel}
 be null for any non null integer value, and must value $1$ for $k=0$ :
 
 \begin{equation}
-    \forall x \in \ZZ :  \KernI(k) =   \delta_0 \label{KernIntDelta0}
+    \forall k \in \ZZ :  \KernI(k) =   \delta_0 \label{KernIntDelta0}
 \end{equation}
 
 For practicall reason, we cannot compute infinite sum like in equation~\ref{KernIntDef}, and we 
@@ -526,6 +526,150 @@ \subsection{Tabulated interpolators}
     % = = = = = = = = = = = = = = = = = = = = = = =
     % = = = = = = = = = = = = = = = = = = = = = = =
 
+\section{Choice of an Interpolator}
+
+    % = = = = = = = = = = = = = = = = = = = = = = =
+
+\subsection{Interpolator for image up-izing}
+
+Interpolator are mainly usefull for image resampling.  When the resampling is done
+to a higher resolution, or a resolution equivalent the choice of an interpolator,
+there is no so much to say, and it is mainly an affair of trade-off between accuracy and efficiency :
+
+\begin{itemize}
+   \item for fast computation and low accuracry, bi-linear can be a good choice;
+         it's not recommanded when derivative are required;
+
+   \item bi-cubic can be a good compromize quality/efficiency, the default recommander
+         parameter being $P=-\frac{1}{2}$, but higher value can be used for image enhancing
+         when used at high zooming;
+
+   \item for "best" ressampling, with minimal aliazing, according to signal theory, the sinus cardinal
+         can be used;
+
+   \item finally, in image correlation the non standard interpolator \ref{MMVIIInterpol}, especially $M_2$,
+         can be interesting compromize.
+
+\end{itemize}
+
+    % = = = = = = = = = = = = = = = = = = = = = = =
+
+\subsection{Interpolator for image downsizing}
+
+   %   -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
+
+\subsubsection{No aliasing property}
+
+In image downsizing, we can still use interpolator and formula~\ref{KernIntDef},
+however we need to take some precaution to  avoid aliazing. We make the analyis
+in $1d$, the generalization to $Nd$ being direct.
+
+When computing an image $I_S$ from an image $I$ with a  
+down sizing of size $S$, with a kernem $\KernI_S$, we write formula~\ref{KernIntDef}  this way :
+
+\begin{equation}
+     I_S(k') =   \sum_{k}  I(k)  \KernI_S(S k'-k)   
+\end{equation}
+
+To avoid aliasing, we must take care that each $v_k$ contributes equally to the final result :
+
+\begin{equation}
+    \sum_{k'}  \KernI_S(S k'-k)  = 1 \label{Interp:DS1}
+\end{equation}
+
+If we consider an initial kernel $\KernI$ and define :
+
+\begin{equation}
+    \KernI_S(x) = \frac{\KernI(\frac{x}{S})}{S}
+\end{equation}
+
+
+Then property~\ref{Interp:DS1} is  satistied if initial kernel $\KernI$ is
+a partition of unity :
+
+\begin{equation}
+    \sum_{k'}  \KernI_S(S k'-k)  =  \sum_{k'}  \KernI( k'-\frac{k}{S})  = 1
+\end{equation}
+
+All interpolator used in \PPP are partition of unity (at least once tabulated).
+For downsizing is seems  \emph{"Natural"} to select an interpolator with only positive
+coefficient, it we want to take as model the physcical sensor.  An also a smooth
+function seems preferable :
+
+\begin{itemize}
+    \item $SinC$ and cubic for  $p \neq 0$ are non postive interpolator ;
+    \item  nearest neighboor is certainly not smooth;
+    \item  linear interpolator, cubic interpolator for $p=0$ are continous interporlator;
+    \item  $M_2$ interpolator can also be used, it's draw back is to be non standard,
+           and maybe to have a slight blurring effect
+\end{itemize}
+
+
+   %   -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
+
+\subsubsection{"No spatial bias" property}
+
+Someway,  cubic interpolator  with $p=0$, may seems the best choice : smooth (differentiable), positive and support
+reduced to $[-1,1]$.  But there  is another criteria that we need to take care.   We want that the
+downsampling does not create any spatial bias. By that we define the centroid $C(I)$ of $I$
+as  :
+
+\begin{equation}
+    C(I)  =  \frac{\sum_k k I(k)} {\sum_k I(k)}
+\end{equation}
+
+Considering that $I$ is already normalised such $\sum_k I(k)=1$, we have
+$C(I) = \sum_k k I(k)$. By no spatial bias, we mean formally :
+
+
+\begin{equation}
+    C(I_S)  =  \frac{C(I)}{S} \label{Interp:Eq:Centroid}
+\end{equation}
+
+Due to linearity, it is sufficient that equation \ref{Interp:Eq:Centroid} is satisfied
+for functions kroneckers's function $\delta^x$. We obviously have $C(\delta^x) = x$.
+We have :
+
+\begin{equation}
+    \delta^x_S(k') = \frac{1}{S} *  \sum_{k}  \delta_x(k)  \KernI_S(S k'-k ) =  \frac{1}{S} \KernI_S(S k'-x)
+\end{equation}
+
+So we have :
+
+\begin{equation}
+    C(\delta^x_S) =  \sum_{k'}  \KernI(k'-\frac{x}{S}) k'  
+\end{equation}
+
+
+So we finnally have the equation for non spatially biased interpolator :
+
+\begin{equation}
+   \forall x     \sum_{k'}  \KernI(k'-\frac{x}{S}) k'  = \frac{x}{S} \label{Interp:NoSpatialBias}
+\end{equation}
+
+
+The function {\tt TestBiasInterpolator} of \PPP make an experimental test on formula
+\ref{Interp:NoSpatialBias}. The results of this test is the following :
+
+\begin{itemize}
+    \item  linear, $SinC$ and   $M_2$  are not spatially biased;
+
+    \item  cubic interpolator for $p=-0.5$ is not spatially biased (this is the interpolator that
+           transformat lines in lines);
+
+    \item  other cubic interpolator are spatially biased .
+\end{itemize}
+
+So finally, the cubic interpolator is not such a goof choice for image down sampling, 
+because the only  non negative  is biased.  This used to be the default choice for
+image down sampling, but it's no longer the case. For now the default is the linear
+interpolator, maybe coud be replaced by $M_2$ ?
+
+
+    % = = = = = = = = = = = = = = = = = = = = = = =
+    % = = = = = = = = = = = = = = = = = = = = = = =
+    % = = = = = = = = = = = = = = = = = = = = = = =
+
 \section{Users's view}
 
 \label{InterpUserView}

MMVII/include/MMVII_Image2D.h

@@ -395,7 +395,7 @@ template <class Type>  class cIm2D
 
        /**  a more "sophisticated" version, adapted to non integer factor and using interpolator adapted to scaling */
        cIm2D<Type>  Scale(tREAL8 aFX,tREAL8 aFY=-1,tREAL8 aSzSinC=-1,tREAL8 DilateKernel=1.0,
-                           const std::vector<std::string> & aVNameKernI = {"MMVII"}
+                           const std::vector<std::string> & aVNameKernI = {"Linear"}  // Kernel for down scale
                          ) const;
        /// Version allowing to fix the interpolator (called by version above)
        cIm2D<Type>  Scale(const cInterpolator1D &,tREAL8 aFX,tREAL8 aFY=-1) const;

MMVII/src/Bench/BenchInterpolators.cpp

@@ -3,8 +3,6 @@
 namespace MMVII
 {
 
-
-
 /**  Basic test on bilinear interpolator :   check that "cLinearInterpolator" works,
  *  "GetValueInterpol" work, and "Tabulated Interpolator" works.
  *

MMVII/src/UtiMaths/Interpolators.cpp

@@ -791,8 +791,87 @@ tREAL8  cMultiScaledInterpolator::DiffWeight(tREAL8  anX) const
      return aRes;
 }
 
+/* **************************************************** */
+/*                                                      */
+/*          Test bias on centroid of interpolator       */
+/*                                                      */
+/* **************************************************** */
+
+tREAL8  CentroidDownScale_Interpolator(const cDiffInterpolator1D & anInt,tREAL8 aScale,tREAL8 aPhase)
+{
+
+   //  Down scale is given by
+   //
+   //     IS(x) =  Sum[d]   I(Sx+d) K(d/S)   {1}
+   //    
+   //  Supose that I  is a delta function at phase  A :
+   //    * in {1} the only non null term  is given for Sx+d=A  <=>  d = A-Sx
+   //    * the 
+   //      V(x) = I(A) K(A/S -x)
+   // Now the centroid, C(A)  of V(x) is given by 
+   //     Sum [x] = x I(A) K(A/S -x)
+
+    tREAL8 aAmpl = anInt.SzKernel() +2;
+    tREAL8 aCenter = aPhase / aScale;
+
+    cWeightAv<tREAL8,tREAL8> aAvg;
+    for (int aX=   round_down(aCenter-aAmpl) ; aX<= round_up(aCenter+aAmpl+1e-5) ; aX++)
+    {
+        tREAL8 aW = anInt.Weight(aCenter-aX);
+        if (aW!=0)
+        {
+           aAvg.Add(aW,aX);
+           // StdOut()  << " X=" << aX  << " C-X=" << aCenter-aX<< " W=" << aW << "\n";
+        }
+    }
+    // StdOut() << "CCIII=============== " << aCenter << "\n";
+
+    return aAvg.Average();
+}
+
+tREAL8 TestBiasInterpolator(const cDiffInterpolator1D & anInt,tREAL8 aScale,tREAL8 aP0)
+{
+    tREAL8 aSumDif=0;
+    for (int aK= -aScale-1; aK<=aScale+1 ; aK++)
+    {
+         tREAL8 aC = CentroidDownScale_Interpolator(anInt,aScale,aK+aP0);
+         aSumDif += std::abs(aC -(aK+aP0)/aScale);
+    }
+    StdOut() << " INTERP BIAS=" << aSumDif << " For Scale=" << aScale << " Phase=" << aP0<< "\n";
+    return aSumDif;
+}
+
+
+void TestBiasInterpolator(const cDiffInterpolator1D & anInt)
+{
+   StdOut() << "=====  INTERP : " << anInt.VNames() << "\n";
+   TestBiasInterpolator(anInt,3.0,0.0);
+   TestBiasInterpolator(anInt,5.0,0.0);
+   TestBiasInterpolator(anInt,5.5,0.0);
+   TestBiasInterpolator(anInt,5.0,0.1);
+   TestBiasInterpolator(anInt,5.0,0.5);
+}
+
+void TestBiasInterpolator()
+{
+    //  Not OK
+    TestBiasInterpolator(cCubicInterpolator(-1.0));
+    TestBiasInterpolator(cCubicInterpolator(0.0));
+
+    //  OK
+    TestBiasInterpolator(cLinearInterpolator());
+    TestBiasInterpolator(cMMVII2Inperpol());
+    TestBiasInterpolator(cSinCApodInterpolator(20,20));
+
+    TestBiasInterpolator(cCubicInterpolator(-0.5));
+}
+
+/* **************************************************** */
+/*                                                      */
+/*          Test scale image vs centroid                */
+/*                                                      */
+/* **************************************************** */
 
-// tREAL8  cMultiScaledInterpolator::DiffWeight(tREAL8  anX) const { return 0.0; }
 
 template <class Type> cPt2dr  Centroid(const cDataIm2D<Type> & aDIm)
 {
@@ -804,33 +883,59 @@ template <class Type> cPt2dr  Centroid(const cDataIm2D<Type> & aDIm)
   return aAvg.Average();
 }
 
-template <class Type> void BenchScaleIm(Type aValC,tREAL8 aEps)
+template <class Type> void BenchScaleIm(Type aValC,tREAL8 aEps,const cDiffInterpolator1D & aInt)
 {
     for (int aK=0 ; aK<50; aK++)
     {
           tREAL8 aScale = 2.0 +(aK%4) ;
 
-          cPt2di aSz(round_ni(RandInInterval(30,40)),round_ni(RandInInterval(30,40)));
+          cPt2di aSz(round_ni(RandInInterval(50,60)),round_ni(RandInInterval(50,60)));
           cPt2di aMil = aSz/2;
 
           cIm2D<Type> aIm(aSz,nullptr,eModeInitImage::eMIA_Null);
           aIm.DIm().SetV(aMil,aValC);
 
-          cIm2D<Type>  aImSc  = aIm.Scale(aScale);
-
+          std::unique_ptr<cTabulatedDiffInterpolator>  aTabInt ( cScaledInterpolator::AllocTab(aInt,aScale,1000));
+          cIm2D<Type>  aImSc  = aIm.Scale(*aTabInt,aScale);
+
 
           tREAL8 aDif = Norm2(Centroid(aIm.DIm())  - Centroid(aImSc.DIm()) * aScale);
-          // StdOut() << " CC  " <<  aDif << "\n";
-          MMVII_INTERNAL_ASSERT_bench(aDif<aEps,"BenchScaleIm");
+          if (aDif>=aEps)
+          {
+              StdOut() << " DIF= " << aDif << " Sc=" << aScale 
+                       << " x:" << aMil.x() %(int)aScale 
+                       << " y:" << aMil.y() %(int)aScale 
+                       << "\n";
+              MMVII_INTERNAL_ASSERT_bench(aDif<aEps,"BenchScaleIm");
+          }
     }
 }
 
+
+
 void Bench_cMultiScaledInterpolator()
 {
+     // TestBiasInterpolator();
+
+     // just to test that we dont have problem with default interpolator
+     {
+         cIm2D<tREAL4> aIm(cPt2di(100,100));
+         aIm.Scale(10);
+     }
 
-     BenchScaleIm<tREAL4>(1.0,1e-5);
-     BenchScaleIm<tU_INT2>(50000,1e-2);
-     BenchScaleIm<tINT4>(100000000,1e-4);
+     if (1)
+     {
+        // Dont understand why bicubic[-0.5]  works except sometimes ...
+        BenchScaleIm<tREAL8>(1.0,1e-5,cCubicInterpolator(-0.5));
+     }
+     BenchScaleIm<tREAL8>(1.0,1e-2,cCubicInterpolator(-0.5));
+
+     BenchScaleIm<tREAL8>(1.0,1e-5,cLinearInterpolator());
+
+     BenchScaleIm<tREAL4>(1.0,1e-5,cMMVII2Inperpol());
+     BenchScaleIm<tU_INT2>(50000,1e-2,cMMVII2Inperpol());
+     BenchScaleIm<tINT4>(100000000,1e-4,cMMVII2Inperpol());
+     // BenchScaleIm<tREAL8>(1.0,1e-5,cSinCApodInterpolator(10,10));
      // BenchScaleIm<tREAL8>(1.0);
 
 
@@ -851,16 +956,6 @@ void Bench_cMultiScaledInterpolator()
               aMSI.SetScale(aS);
           }
      }
-
-     /*
-     cPt2di  aSz(100,120);
-     cPt2di  aPMil = aSz/2;
-     cIm2d<tU_INT1>  aImDirac(aSz);
-     aImDirac.aImDirac(aPMil,1);
-     */
 }
 
-
-
-
 };