--- a
+++ b/Thoracic Organs Segmentation code/OpticalFlow3D/OpticalFlow3D.cpp
@@ -0,0 +1,158 @@
+/*-----------------------------------------------------------------------
+# This is a modified version of the following file:
+#  
+#   File        : optical_flow.cpp 
+#   Description : Compute the optical flow between two images, with a multiscale and variational algorithm  
+#   Author      : Tschumperlé David 
+#   Institution : ODYSSEE, INRIA Sophia Antipolis. 
+#   Contact     : David.Tschumperle@sophia.inria.fr 
+#   Date        : 03/12/2003 
+#    
+#   This program is free software; you can redistribute it and/or modify 
+#   it under the terms of the GNU General Public License as published by 
+#   the Free Software Foundation; either version 2 of the License, or 
+#   (at your option) any later version. 
+//////////////////////////////////////////////////////////////////////////////////////////
+#The Input parameters are two images (2D or 3D) and as an output we have the displacement vector u.
+--------------------------------------------------------------------------*/
+
+#include <mex.h>
+#include <mat.h>
+#include <matrix.h>
+
+#define cimg_plugin "cimgmatlab.h"
+
+#include "CImg.h"
+#include <iostream>
+#include <string>
+#include <fstream>
+
+using namespace cimg_library;
+using namespace std;
+
+
+//global constants
+const float smooth = 0.1f;//"Flow Smoothness"
+const float precision = 0.09f;//"Convergence precision"
+const unsigned int nb=110;//"Number of warped frames"
+const unsigned int nbscale = 0 ;//"Number of scales (0=auto)");
+
+const bool normalize = true; //"Histogram normalization of the images"
+const bool morph = true;//"Morphing mode"
+const bool imode  = true;//"Complete interpolation (or last frame is missing)"
+const bool dispflag = true;//"Visualization"
+//Functions
+CImg<> optmonoflow(const CImg<>& I1, const CImg<>& I2, const CImg<>& u0,
+                   const float smooth, const float precision, CImgDisplay& disp);
+CImg<> optflow(const CImg<>& xsrc, const CImg<>& xdest,
+               const float smooth, const float precision, const unsigned int pnb_scale, CImgDisplay& disp);
+                
+//mex function
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+  if (nrhs < 2) mexErrMsgTxt("No enough input arguments.");
+  if (nrhs > 2) mexErrMsgTxt("Too many input arguments.");
+  if (nrhs == 2){
+                CImg<> src_blur(prhs[0],false), dest_blur(prhs[1],false);//2 Input Volumes
+                //Input images preprocessing
+                // src_blur = normalize?src_blur.get_blur(0.5f).equalize(256):src_blur.get_blur(0.5f),
+               //  dest_blur  = normalize?dest_blur.get_blur(0.5f).equalize(256):dest_blur.get_blur(0.5f);
+                 CImgDisplay disp;
+                 const CImg<> u = optflow(src_blur,dest_blur,smooth,precision,nbscale,disp);
+                
+                 plhs[0] = u.toMatlab();
+                 ofstream outputFilex, outputFiley,outputFilez;
+                 outputFilex.open("displacement_x.txt");
+                 outputFiley.open("displacement_y.txt");
+                 outputFilez.open("displacement_z.txt");
+                 cimg_forY(u,y){
+                      cimg_forXZ(u,x,z){
+                    //float mag = (float) sqrt(pow(u(x,y,z,0),2)+pow(u(x,y,z,1),2)+pow(u(x,y,z,2),2));
+                          outputFilex<<u(x,y,z,0)<<" ";
+                          outputFiley<<u(x,y,z,1)<<" ";
+                          outputFilez<<u(x,y,z,2)<<" ";
+                     }
+                  outputFilex<<endl;outputFiley<<endl;outputFilez<<endl;
+                 }
+                outputFilex.close();outputFiley.close();outputFilez.close();
+            
+ }
+  return;
+     
+}
+/////////////////////////////////////////////////////////////////////////////////
+
+CImg<> optflow(const CImg<>& xsrc, const CImg<>& xdest,
+               const float smooth, const float precision, const unsigned int pnb_scale, CImgDisplay& disp) {
+  const CImg<>
+    src  = xsrc.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),xdest.dimz()).normalize(0,1),
+    dest = xdest.get_pointwise_norm(1).resize(xdest.dimx(),xdest.dimy(),xdest.dimz()).normalize(0,1);
+    CImg<> u =  CImg<>(src.dimx(),src.dimy(),src.dimz(),3,1).fill(0);
+       
+  //multi-scalling 
+    const unsigned int nb_scale = pnb_scale>0?pnb_scale:(unsigned int)(2*log((double)(max(max(src.dimx(),src.dimy()),src.dimz()))));
+  //const unsigned int nb_scale = pnb_scale>0?pnb_scale:(unsigned int)(2*std::log((double)(cimg::max(src.dimx(),src.dimy()))));
+  for (int scale=nb_scale-1; scale>=0; scale--) {
+       const CImg<> I1 = src.get_resize((int)(ceil(src.dimx()/std::pow(1.5,scale))), (int)(ceil(src.dimy()/std::pow(1.5,scale))) ,(int)(ceil(src.dimz()/std::pow(1.5,scale))));
+       const CImg<> I2 = dest.get_resize((int)(ceil(src.dimx()/std::pow(1.5,scale))), (int)(ceil(src.dimy()/std::pow(1.5,scale))) ,(int)(ceil(src.dimz()/std::pow(1.5,scale))));
+       u*=1.5;
+       u = optmonoflow(I1,I2,u,smooth,(float)(precision/std::pow(2.25,1+scale)),disp);//apply optical flow algorithm for different scales
+    
+  }
+    cimg_forXYZV(u,x,y,z,v){
+    if(src(x,y,z)==0) u(x,y,z,v)=0;}
+ return u;
+}
+
+// optmonoflow() : estimate optical flow for one scale ( semi-implicite PDE scheme ) between I2->I1
+//---------------
+CImg<> optmonoflow(const CImg<>& I1, const CImg<>& I2, const CImg<>& u0,
+                   const float smooth, const float precision, CImgDisplay& disp) {
+
+ CImg<float> u = u0.get_resize(I1.dimx(),I1.dimy(),I1.dimz(),3);
+ CImg<float> dI(I2.dimx(),I2.dimy(),I2.dimz(),3);
+  dI.fill(0);
+ //CImg_3x3x3(I,float);
+  float dt=2,E=1e20f;
+
+  // compute first derivatives of I2
+    cimg_for3XYZ(dI,x,y,z){
+              dI(x,y,z,0) = 0.5*(I2(_n1x,y,z)-I2(_p1x,y,z));//derivatives of I2
+              dI(x,y,z,1) = 0.5*(I2(x,_n1y,z)-I2(x,_p1y,z));
+              dI(x,y,z,2) = 0.5*(I2(x,y,_n1z)-I2(x,y,_p1z));
+                          
+     }
+   
+// Main PDE iteration
+for (unsigned int iter=0; iter<50; iter++) {
+    std::fprintf(stderr,"\r- Iteration %d - E = %g",iter,E); std::fflush(stderr);
+    const float Eold = E;
+    E = 0;
+    cimg_for3XYZ(u,x,y,z) { //3x3 neighborhood  
+      const float
+        X = x + u(x,y,z,0),
+        Y = y + u(x,y,z,1),
+        Z = z + u(x,y,z,2),
+        deltaI = (float)(I2.linear_atXYZ(X,Y,Z) - I1(x,y,z));// partial derivative over time-dI/dt
+      float tmpf = 0;
+      cimg_forV(u,k) {
+        const float
+          ux  = 0.5f*(u(_n1x,y,z,k)-u(_p1x,y,z,k)),//derivatives of velocity(displacemenet)
+          uy  = 0.5f*(u(x,_n1y,z,k)-u(x,_p1y,z,k)),
+          uz  = 0.5f*(u(x,y,_n1z,k)-u(x,y,_p1z,k)); 
+          u(x,y,z,k) = (float)( u(x,y,z,k) +
+                            dt*(
+                                -deltaI*dI.linear_atXYZ(X,Y,Z,k) +
+                                smooth* ( u(_n1x,y,z,k) + u(_p1x,y,z,k) + u(x,_n1y,z,k) + u(x,_p1y,z,k)+  u(x,y,_n1z,k) + u(x,y,_p1z,k)-6*u(x,y,z,k))
+                                )/(1+4*smooth*dt)
+                              );
+      tmpf += ux*ux + uy*uy + uz*uz;
+     
+      }
+      E += deltaI*deltaI + smooth * tmpf;
+    }
+    if (cimg::abs(Eold-E)<precision) break;
+    if (Eold<E) dt*=0.5;
+ 
+}
+  return  u;
+}
\ No newline at end of file