Transformation.cpp

Go to the documentation of this file.

/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * 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.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include<ImLib3D/FFT.hpp>
#include<ImLib3D/Transformation.hpp>
#include<ImLib3D/Threshold.hpp>
#include<ImLib3D/Arithmetic.hpp>
#include<ImLib3D/Convolution.hpp>



template <>
void
IP3D::TransformWithInverseField<Mask3D>(const Mask3D& src, const Field3Df& field, Mask3D& res0)
{
    // Image Args
    ImageProcessorArgs<Mask3D> args(ImArgs(src),res0);
    args.SetDestSize(field.Size());
    Mask3D &res=args.GetDest();

    // find distinct values
    vector<int> values0;
    values0.assign(256,0);
    for(Mask3D::const_iteratorFast mp=src.begin();mp!=src.end();++mp)
    {
        values0[*mp]++;
    }
    vector<byte> distinct;
    uint i;
    for(i=0;i<256;i++){if(values0[i]){distinct.push_back(i);}}

    Image3Df resWeight(res.Size());
    resWeight.Fill(.01);
    res.Fill(0);
    // extract mask for each distinct  value
    for(i=0;i<distinct.size();i++)
    {
        byte value=distinct[i];
        Mask3D srcMask;
        IP3D::SimpleThresholds2(src,value,value,srcMask);
        srcMask.Normalize();
        Image3Df ima;
        // transform mask with field and good interpol
        {
            ima=srcMask;
            // smooth mask to avoid violent aliasing
            Filter filter(3,3,3,Vect3Di(1,1,1));
            filter.Fill(1.0/filter.GetNVoxels());
            IP3D::Convolution(ima,filter,ima);
            // apply transform 
            IP3D::TransformWithInverseField(ima,field,ima);
        }
        // add result to final mask
        {
            Mask3D::iteratorFast rp;
            Image3Df::iteratorFast ip;
            Image3Df::iteratorFast wp;
            // use weight image to vote for best final pixel value
            for(rp=res.begin(),ip=ima.begin(),wp=resWeight.begin();
                rp!=res.end();
                rp++,ip++,wp++)
            {
                if(*ip>*wp)
                {
                    *rp=value;
                    *wp=*ip;
                }
            }
        }
//          mprintf("** affine trf image:for value:%d\n",value);ImLib3DDisplay::Show(ima);
//          mprintf("** weights after :%d\n",value);ImLib3DDisplay::Show(resWeight);
//          mprintf("** res after :%d\n",value);ImLib3DDisplay::Show(res);
        

    }
}


//  void 
//  CTransform::CreateField(Field3Df &field, const Affine3DTransform& transfo)
//  {
//      Matrix M = transfo.M;
//      ColumnVector T = transfo.T;
//      float a11=M(1,1),a21=M(2,1),a31=M(3,1);
//      float a12=M(1,2),a22=M(2,2),a32=M(3,2);
//      float a13=M(1,3),a23=M(2,3),a33=M(3,3);

//      Field3Df::iteratorFast p;
//      p=field.begin();
//      Vect3Df T1;
//      Vect3Df Ttemp(T(1), T(2), T(3));
//      for(int z=0;z<field.Depth() ;z++)
//      {
//          Vect3Df Vz(a13,a23,a33);
//          Vz*=z;
//          for(int y=0;y<field.Height();y++)
//          {
//              Vect3Df Vy(a12,a22,a32);
//              Vy*=y;
//              Vy+=Vz;
//              for(int x=0;x<field.Width() ;x++)
//              {
//                  Vect3Df Vx(a11,a21,a31);
//                  Vx*=x;
//                  *p=Vy+Vx+Ttemp-Vect3Df(x,y,z);
//                  p++;
//              }
//          }
//      }
//  }

void 
IP3D::TransformAffineHelpers::CreateInverseField(Field3Df &field, const Affine3DTransform& transfo0)
{
//      cout << "CreateField:" << endl << transfo0.M << "T:"<< transfo0.T << endl;
    Affine3DTransform transfo=transfo0;
    transfo.Inverse();
//      cout << "CreateField:" << endl << transfo.M << "T:"<< transfo.T << endl;
    Field3Df::iteratorXYZ fp;
    ImageProgress tracker("CTransformAffine::CreateInverseField",fp);
    for(fp=field.begin();fp!=field.end();++fp)
    {
        Vect3Df P(fp.x,fp.y,fp.z);
        *fp=transfo(P)-P;
    }
//      Matrix M = transfo.M;
//      ColumnVector T = transfo.T;
//      float a11=M(1,1),a21=M(2,1),a31=M(3,1);
//      float a12=M(1,2),a22=M(2,2),a32=M(3,2);
//      float a13=M(1,3),a23=M(2,3),a33=M(3,3);

//      Field3Df::iteratorFast p;
//      p=field.begin();
//      Vect3Df T1;
//      Vect3Df Ttemp(T(1), T(2), T(3));
//      for(int z=0;z<field.Depth() ;z++)
//      {
//          Vect3Df Vz(a13,a23,a33);
//          Vz*=z;
//          for(int y=0;y<field.Height();y++)
//          {
//              Vect3Df Vy(a12,a22,a32);
//              Vy*=y;
//              Vy+=Vz;
//              for(int x=0;x<field.Width() ;x++)
//              {
//                  Vect3Df Vx(a11,a21,a31);
//                  Vx*=x;
//                  *p=Vy+Vx+Ttemp-Vect3Df(x,y,z);
//                  p++;
//              }
//          }
//      }

}

//  template <>
//  void IP3D::TransformWithFieldInt<Mask3D>(const Mask3D& imOrig, const Field3Df& field, Mask3D& _imRes)
//  {
//      Image3Df src1;src1=imOrig;
//      Image3Df res1;
//      Apply( src1, field, res1);
//      res1+=.5;
//      float minv,maxv;
//      res1.FindMinMax(minv,maxv);mprintf("** minv:%f,maxv:%f\n",minv,maxv);
//      IP3D::LimitThreshold(res1,0.0f,255.0f,res1);
//      res1.FindMinMax(minv,maxv);mprintf("** minv:%f,maxv:%f\n",minv,maxv);
//      _imRes=res1;
//  }

---