/tomo/pyhst

#include"Python.h"
#include"structmember.h"
#include<string.h>

#include "numpy/arrayobject.h"

#ifdef LINUX
 #define _FILE_OFFSET_BITS  64
#endif


#define DEBUG(a)


 /* The error object to expose 
 */

static PyObject *ErrorObject;
#define onError(message)\
  { PyErr_SetString(ErrorObject, message); return NULL;}


#define  __CCD_Filter_name__  "CCD_FilterAbsRel"
#define __init__Filter__      initCCD_FilterAbsRel

static char Filter_doc[]=
"/* "
""
"   Filter( item , pos_edf_, size_edf_, pos_edf, size_edf,   CCD_FILTER_PARA     )"
"  "
"   CCD_FILTER_PARA  is a python dictionary with parameters"
" */"
;

static PyObject *
Filter(PyObject *self_a, PyObject *args)
{
  PyObject * A         ;
  PyObject * pos_      ;
  PyObject * siz_      ;
  PyObject * pos       ;
  PyObject * siz       ;
  PyObject * para_dict ;
  PyObject * tmp;
  char msg[400];
  
  int x0_,x1_, x0,x1, s0_,s1_,s0,s1;
  int X0,X1,X0_, X1_;
  int dim0, dim1;
  int nsum;
  int ds;


  int c0,c1;
  int c0_,c1_;
  float temp[9];
  float swap;
  int bool;
  int count;


  PyObject * dic_key;
  PyObject * dic_value;

  double threshold ;
  int fracornot;
  int absornot;

  double valtocomp;
 

  int i;


  int *iptr[] = { &x0_ , &x1_  ,&x0 , &x1 , &s0_ , &s1_ , &s0 , &s1 };

  int shift_x[] = { -1,-1,-1, 0, 0,  0, 1,  1, 1  };
  int shift_y[] = { -1, 0, 1,-1, 0, 1,-1,  0, 1  };
  int si;

  PyObject ** containers[] = {&pos_ , &pos, &siz_, &siz   };


  float * a_data;
  float * m_data;





 if(!PyArg_ParseTuple(args,"OOOOOO:CDD_Filter",&A,&pos_,&siz_, &pos, &siz, &para_dict  ))
   return NULL;

 /* check the Objects */

  if(!PyArray_Check(A ))    onError("not a PyArray, argument 1");
  if(!PyList_Check (pos_ ) )     onError("not a PyList, argument 2");
  if(!PyList_Check (siz_ ))    onError("not a PyList, argument 3");
  if(!PyList_Check (pos ))         onError("not a PyList, argument 4");
  if(!PyList_Check(siz ))       onError("not a PyList, argument 5");
  if(!PyDict_Check(para_dict ))        onError("not a PyDict, argument 6");

  /* check the types */

  if( ((PyArrayObject *) A )  ->descr->type_num != PyArray_FLOAT ) onError(" arg 1 is not an array of float " ) ;

  /* check n of dimensione */

  if( ((PyArrayObject *) A)->nd !=2)  onError(" arg 1 is not a matrix " ) ;
 

  /* check contiguity */

  if(  ((PyArrayObject *) A) ->flags %2 == 0) onError(" array has to be contiguous");

  a_data = (float *) ((PyArrayObject *) A)->data;
  /* getting positions and sizes */

  for (i=0; i<8; i++) {
    tmp =  PyList_GetItem( *(containers[i/2])  , i%2);
    if(!PyInt_Check( tmp )) {
      sprintf(msg, " not an integer, item %d of argument %d", i%2, i/2+2);
      onError(msg);
    }
    *iptr[i] = PyInt_AS_LONG( tmp ); 
  }


  /* getting parameter(s) contained in the dictionary */

  dic_key = PyString_FromString("threshold");
  dic_value =  PyDict_GetItem(para_dict , dic_key);
  Py_DECREF(dic_key);

  if( dic_value == NULL) 
    onError(" unable to find the threshold key in the passed dictionary");

  if( !PyFloat_Check(dic_value) ) 
    onError("supplied value for threshold is not a float ");
  
  threshold = PyFloat_AsDouble (dic_value  );



  /* ABS???????????????????/ */
  dic_key = PyString_FromString("ABS");
  dic_value =  PyDict_GetItem(para_dict , dic_key);
  Py_DECREF(dic_key);
  if( dic_value == NULL) 
    onError(" unable to find the ABS key in the passed dictionary");
  if( !PyInt_Check(dic_value) ) 
    onError("supplied value for ABS  is not a integer ");
  
  absornot = PyInt_AsLong  (dic_value  );
  if( absornot !=0 && absornot!=1 ) {
    onError("supplied value for ABS  must be 1 or 0 ");
  }
  
  
  
  /* FRACTIONNAIRE???????????????????/ */
  dic_key = PyString_FromString("FRAC");
  dic_value =  PyDict_GetItem(para_dict , dic_key);
  Py_DECREF(dic_key);
  if( dic_value == NULL) 
    onError(" unable to find the FRAC key in the passed dictionary");
  if( !PyInt_Check(dic_value) ) 
    onError("supplied value for  FRAC is not a integer ");
  
  fracornot = PyInt_AsLong  (dic_value  );
  if( fracornot !=0 && fracornot !=1 ) {
    onError("supplied value for  FRAC  must be 1 or 0 ");
  }
  
  
  /* setting better the limits */

  x0 = x0 - x0_ ;
  x1 = x1 - x1_ ;
  
  x0_ = 0 ;
  x1_ = 0 ;
  

  X0 = x0 + s0 ;
  X1 = x1 + s1 ;
  
  X0_ = x0_ + s0_ ;
  X1_ = x1_ + s1_ ;


  /* checking that the supllied limits are not completely wrong */

  if ( ((PyArrayObject *) A)->dimensions[0] != s0_  || ((PyArrayObject *) A)->dimensions[1] != s1_  ) {
    onError(" array dimension and supplied limits are not consistent");
  }
  
  dim0 = s0_;
  dim1 = s1_;
  
  /* straighforward addressing scheme in case of contiguous memory. Or you'll need to use strides */


  m_data=(float*) malloc( dim0*dim1 * sizeof(float) ) ;
  memcpy(m_data ,a_data  , dim0*dim1*sizeof(float) );
  
  
#define AA(i,j)  a_data[(i)*dim1+(j)    ]
#define MM(i,j)  m_data[(i)*dim1+(j)    ]
  
  for(c0=x0+1 ; c0< X0 -1; c0++) {
    
    for(c1=x1+1; c1< X1 -1; c1++) {
      
      
      //Naechste Nachbarn einlesen
      temp[0] = AA(c0-1,c1-1);
      temp[1] = AA(c0-1,c1);
      temp[2] = AA(c0-1,c1+1);
      temp[3] = AA(c0,c1-1);
      temp[4] = AA(c0,c1);
      temp[5] = AA(c0,c1+1);
      temp[6] = AA(c0+1,c1-1);
      temp[7] = AA(c0+1,c1);
      temp[8] = AA(c0+1,c1+1);
      
      //Bubble-Sort
      do
	{
	  bool = 0;
	  for (count = 0; count < 8; count++)
	    {
	      if (temp[count] > temp[count+1])
		{
		  swap = temp[count];
		  temp[count] = temp[count+1];
		  temp[count+1] = swap;
		  bool = 1;
		}
	    }
	} while(bool!=0);

      valtocomp =       AA(c0,c1) -temp[4] ;
      if( absornot){
	valtocomp = fabs( valtocomp);
      }
      if( fracornot ) {
	valtocomp /= fabs(AA(c0,c1));
      }
      if (  valtocomp >= threshold ) MM(c0,c1) = temp[4];

      else MM(c0,c1) = AA(c0,c1);
    }
  }
  
  ds=(s0-1);
  if(ds==0) ds=1;
  for(c0= x0 ; c0<X0 ; c0 +=ds) {  
    for(c1=x1; c1<X1; c1 ++ ) {   
      nsum=0;
      
      for(si=0; si<9; si++) {
	
	c0_ = c0+ shift_x[si];
	c1_ = c1+ shift_y[si];
	
	
	
	if( c0_>=x0_ &&  c0_ <s0_ &&    c1_>=x1_ &&  c1_ <X1_  ) {
	  temp[nsum ] = AA(c0_,c1_); ;
	  nsum++;
	}
      }
      
      //Bubble-Sort
      do
	{
	  bool = 0;
	  for (count = 0; count < nsum-1; count++)
	    {
	      if (temp[count] > temp[count+1])
		{
		  swap = temp[count];
		  temp[count] = temp[count+1];
		  temp[count+1] = swap;
		  bool = 1;
		}
	    }
	} while(bool!=0);
      
      if (  AA(c0,c1) -temp[nsum/2] >= threshold ) MM(c0,c1) = temp[nsum/2];
      else MM(c0,c1) = AA(c0,c1);
      
    }
    
  }
  for(c0= x0 ; c0<X0 ; c0 ++) {
    ds=(s1-1);
    if(ds==0) ds=1;

    for(c1=x1; c1<X1; c1 += ds) {
      nsum=0;

      for(si=0; si<9; si++) {

        c0_ = c0+ shift_x[si];
        c1_ = c1+ shift_y[si];



        if( c0_>=x0_ &&  c0_ <s0_ &&    c1_>=x1_ &&  c1_ <X1_  ) {
          temp[nsum ] = AA(c0_,c1_); ;
          nsum++;
        }
      }

      //Bubble-Sort
      do
        {
          bool = 0;
          for (count = 0; count < nsum-1; count++)
            {
              if (temp[count] > temp[count+1])
                {
                  swap = temp[count];
                  temp[count] = temp[count+1];
                  temp[count+1] = swap;
                  bool = 1;
                }
            }
        } while(bool!=0);

      if (  AA(c0,c1) -temp[nsum/2] >= threshold ) MM(c0,c1) = temp[nsum/2];
      else MM(c0,c1) = AA(c0,c1);

    }

    
    
  } 
  
  
  memcpy( a_data, m_data, dim0*dim1*sizeof(float) );
  free(m_data);
  Py_INCREF(Py_None);
  return Py_None;
    
}



static PyMethodDef CCD_Filter_functions [] = {
  {"Filter", Filter,  METH_VARARGS, Filter_doc},
  { NULL, NULL}
};






void 
__init__Filter__   ()
{
  PyObject *m, *d;
  m = Py_InitModule(  __CCD_Filter_name__    ,  CCD_Filter_functions  );
  d = PyModule_GetDict(m);
  ErrorObject = Py_BuildValue("s",  __CCD_Filter_name__  ".error");
  PyDict_SetItemString(d,"error", ErrorObject);
  if(PyErr_Occurred())
    Py_FatalError("can't initialize module " __CCD_Filter_name__  );

#ifdef import_array
  import_array();
#endif
}