Commits

Anonymous committed e02a5dc

changes to fluid kernel when gamma is zero

  • Participants
  • Parent commits a967dc3

Comments (0)

Files changed (5)

File Code/Cxx/inc/alg/FluidKernelFFTBase.h

 
 #define MAX_FFT_TABLE_LENGTH 512
 
+// if REPLACE_ZERO_FREQ is defined, FluidKernelFFT will replace the
+// zero frequency when gamma is zero (ie L is non-invertible),
+// otherwise the zero frequency is set to zero
+// #define REPLACE_ZERO_FREQ
+
 namespace PyCA {
 
 // forward declaration

File Code/Cxx/src/alg/CFluidKernelFFT.cxx

     T L10, L11;
     T L20, L21, L22;
 
+#if defined(REPLACE_ZERO_FREQ)
     ComplexT<T> zeroFreqX = this->mFFTArrayX[0];
     ComplexT<T> zeroFreqY = this->mFFTArrayY[0];
     ComplexT<T> zeroFreqZ = this->mFFTArrayZ[0];
+#else
+    ComplexT<T> zeroFreqX = {static_cast<T>(0),static_cast<T>(0)};
+    ComplexT<T> zeroFreqY = {static_cast<T>(0),static_cast<T>(0)};
+    ComplexT<T> zeroFreqZ = {static_cast<T>(0),static_cast<T>(0)};
+#endif
 
     for (int z = 0; z < this->mComplexSize.z; ++z)
 	{

File Code/Cxx/src/alg/GFluidKernelFFT.cu

     // element of FFT arrays) and replace it after kernel
     ComplexT<T> zeroFreqX={0.f,0.f}, 
        zeroFreqY={0.f,0.f}, zeroFreqZ={0.f,0.f};
+#if defined(REPLACE_ZERO_FREQ)
     if (inverseOp && this->mGamma == static_cast<T>(0)){
        cudaMemcpy(&zeroFreqX,this->mFFTArrayX,
 		  sizeof(ComplexT<T>),cudaMemcpyDeviceToHost);
        cudaMemcpy(&zeroFreqZ,this->mFFTArrayZ,
 		  sizeof(ComplexT<T>),cudaMemcpyDeviceToHost);
     }
+#endif
 
     if (inverseOp == true) {
         if (this->mDivergenceFree == true) {

File Code/Python/Display/PyCADisplay.py

         
 def GridPlot(vf, sliceIdx=None, dim='z', every=1,
              isVF=True, color='g', plotBase=True, colorbase='#A0A0FF'):
-    sliceArr = common.ExtractSliceArrVF(vf, sliceIdx, dim)
+    sliceArr = np.squeeze(common.ExtractSliceArrVF(vf, sliceIdx, dim))
     sz = sliceArr.shape
     hID = np.mgrid[0:sz[0], 0:sz[1]]
     d1 = np.squeeze(hID[1, ::every, ::every])
         # force redraw
         plt.draw()
 
+def JacDetPlot(vf, title='Jac. Det.',
+               sliceIdx=None, dim='z',
+               isVF=True, cmap='jet',
+               newFig=False):
+    h = core.ManagedField3D(h.grid(), h.memType())
+    jacdet = core.ManagedImage3D(h.grid(), h.memType())
+    core.Copy(h, vf)
+    if isVF:
+        core.VtoH(h)
+    core.JacDetH(jacdet, h)
+    DispImage(jacdet, title=title,
+              sliceIdx=sliceIdx, dim=dim,
+              cmap=cmap, newFig=newFig)
         
 def SaveSlice(fname, im, sliceIdx=None, dim='z',
               cmap='gray', rng=None, t=True):

File Examples/ElastReg.py

         diffOp = FluidKernelFFTGPU()
 
     # initialize some vars
-    zerovec = Vec3Df(0.0, 0.0, 0.0)
-
     nScales = len(scales)
     scaleManager = MultiscaleManager(origGrid)
     for s in scales:
         # initialize / upsample deformation
         if scaleManager.isFirstScale():
             u.setGrid(curGrid)
-            SetMem(u,zerovec)
+            SetMem(u, 0.0)
         else:
             resampler.updateVField(u)
 
             Mul_I(scratchV, diff)
 
             diffOp.applyInverseOperator(gU, scratchV)
+            # why is this negative necessary?
             MulC_I(gU, -1.0)
 
             # u =  u*(1-VFC*ustep) + (-2.0*ustep)*gU
                     scales = [2,1], 
                     nIters = [1000,1000], 
                     ustep = [0.1, 0.1], 
-                    # fluidParams = [0.5, 0.5, 0.001], 
-                    fluidParams = [0.5, 0.5, 0.0], 
+                    # fluidParams = [1.0, 1.0, 0.001], 
+                    fluidParams = [0.5, 0.0, 0.001], 
+                    # fluidParams = [0.5, 0.5, 0.0], 
                     VFC = 0.2, 
                     plotEvery = 500)