filter_morph.h

/***************************************************************************
                          filter_morph.h  -  description
                             -------------------
    begin                : Fr Nov 29 2019
    copyright            : (C) 2019-2021 by Thies Jochimsen
    email                : thies@jochimsen.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef FILTER_MORPH_H
#define FILTER_MORPH_H
 
#include <odindata/filter_step.h>
#include <odindata/utils.h> // for on_grid(...)
 
 
enum morphOp { erode=0, dilate };
 
 
template <morphOp morphop>
class FilterMorph: public FilterStep {
 
  LDRfloat radius;
 
  STD_string label() const {if(morphop==erode) return "erode"; else return "dilate";}
  STD_string description() const {return label()+" image using spherical kernel as structuring element";}
  bool process(Data<float,4>& data, Protocol& prot) const {
 
    Log<Filter> odinlog(c_label(),"process");
 
    Range all=Range::all();
 
    TinyVector<int,4> shape=data.shape();
    TinyVector<int,3> spatshape(shape(sliceDim), shape(phaseDim), shape(readDim));
 
 
 
    TinyVector<float,3> voxel_spacing;
    voxel_spacing(0)=FileFormat::voxel_extent(prot.geometry, sliceDirection, data.extent(sliceDim));
    voxel_spacing(1)=FileFormat::voxel_extent(prot.geometry, phaseDirection, data.extent(phaseDim));
    voxel_spacing(2)=FileFormat::voxel_extent(prot.geometry, readDirection,  data.extent(readDim));
    ODINLOG(odinlog,normalDebug) << "voxel_spacing=" << voxel_spacing << STD_endl;
 
 
    TinyVector<int,3> enclbox;
    enclbox=radius/voxel_spacing+1;
 
    ODINLOG(odinlog,normalDebug) << "enclbox=" << enclbox << STD_endl;
 
 
    STD_vector<TinyVector<int,3> > neighboffset;
 
    for(int k=-enclbox(0); k<=enclbox(0); k++) {
      for(int j=-enclbox(1); j<=enclbox(1); j++) {
        for(int i=-enclbox(2); i<=enclbox(2); i++) {
          TinyVector<float,3> distvec(voxel_spacing(0)*k, voxel_spacing(1)*j, voxel_spacing(2)*i);
 
          float dist=sqrt(sum(distvec*distvec));
          if(dist<=radius) {
            neighboffset.push_back(TinyVector<int,3>( k, j, i));
          }
        }
      }
    }
 
    unsigned int numof_neigb=neighboffset.size();
    ODINLOG(odinlog,normalDebug) << "numof_neigb=" << numof_neigb << STD_endl;
    for(unsigned int i=0; i<numof_neigb; i++) {
      ODINLOG(odinlog,normalDebug) << "neighboffset[" << i << "]=" << neighboffset[i] << STD_endl;
    }
 
 
    Data<float,3> involdata(spatshape);
    Data<unsigned int,3> neighbcount(spatshape); neighbcount=0;
    Data<float,3> outmask(spatshape); outmask=0.0;
    TinyVector<int,3> index, neighbindex;
 
    for(int irep=0; irep<shape(timeDim); irep++) {
 
      involdata(all,all,all)=data(irep,all,all,all);
 
      for(unsigned int i=0; i<involdata.size(); i++) {
        index=involdata.create_index(i);
 
        if(involdata(index)!=0.0) {
          for(unsigned int j=0; j<numof_neigb; j++) {
            neighbindex=index+neighboffset[j];
 
            if(on_grid<3>(spatshape, neighbindex)) {
              neighbcount(neighbindex)++;
            }
          }
        }
      }
 
      for(unsigned int i=0; i<involdata.size(); i++) {
        index=involdata.create_index(i);
        unsigned int count=neighbcount(index);
         if(morphop==erode && count==numof_neigb) outmask(index)=1.0;
         if(morphop==dilate && count>0) outmask(index)=1.0;
      }
 
      data(irep,all,all,all)=outmask(all,all,all);
 
    }
 
 
    return true;
  }
  FilterStep*  allocate() const {return new FilterMorph<morphop>();}
  void init() {
    radius.set_unit(ODIN_SPAT_UNIT).set_description("radius of kernel");
    append_arg(radius,"radius");
  }
 
};
 
 
 
#endif