From 7bdbd14a14280f774a9f0cbbffc61d71782f9ad0 Mon Sep 17 00:00:00 2001 From: Anna Bergamaschi Date: Wed, 24 May 2023 17:34:02 +0200 Subject: [PATCH] Added the clusterize function for calculating the 2x2 and 3x3 clusters --- src/cluster_reader.c | 38 +++++++++++++++ src/cluster_reader.h | 1 + src/creader_module.c | 107 +++++++++++++++++++++++++++++++++++++++++-- src/data_types.h | 30 ++++++++++++ test.py | 10 ++-- 5 files changed, 178 insertions(+), 8 deletions(-) diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 0d737fd..b71034d 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -43,3 +43,41 @@ int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left){ return nph_read; } +int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout){ + + int32_t tot2[4], t2max; + char c; + int32_t val,tot3; + for (int ic=0; icdata[iy*3+ix]; + tot3+=val; + if (ix<=1 && iy<=1) tot2[0]+=val; + if (ix>=1 && iy<=1) tot2[1]+=val; + if (ix<=1 && iy>=1) tot2[2]+=val; + if (ix>=1 && iy>=1) tot2[3]+=val; + } + } + t2max=tot2[0]; + c=cBottomLeft; + for (int i=1; i<4; i++) { + if (tot2[i]>t2max) { + t2max=tot2[i]; + c=i; + } + } + (cout+ic)->c=c; + (cout+ic)->tot2=t2max; + (cout+ic)->tot3=tot3; + //printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); + + } + return n_clusters; + + +} diff --git a/src/cluster_reader.h b/src/cluster_reader.h index 2ae79d9..827f6d4 100644 --- a/src/cluster_reader.h +++ b/src/cluster_reader.h @@ -6,3 +6,4 @@ //Pure C implementation to read a cluster file int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left); +int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout); diff --git a/src/creader_module.c b/src/creader_module.c index f5477e0..b5445b3 100644 --- a/src/creader_module.c +++ b/src/creader_module.c @@ -28,6 +28,19 @@ static PyArray_Descr *cluster_dt() { dtype_dict = Py_BuildValue( "[(s, s),(s, s),(s, s, (i))]", "x", "u2", "y", "u2","data","i4",9 ); + + PyArray_DescrConverter(dtype_dict, &dtype); + Py_DECREF(dtype_dict); + return dtype; +} + +static PyArray_Descr *cluster_analysis_dt() { + PyObject *dtype_dict; + PyArray_Descr *dtype; + dtype_dict = Py_BuildValue( + "[(s, s),(s, s),(s, s)]", "tot3", "i4","tot2","i4", + "corner", "u4"); + PyArray_DescrConverter(dtype_dict, &dtype); Py_DECREF(dtype_dict); return dtype; @@ -127,11 +140,95 @@ ClusterFileReader_read(ClusterFileReader *self, PyObject *args, return clusters; } + + + + +// read method +static PyObject * +ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args, + PyObject *kwds) +{ + + + + + //Create an uninitialized numpy array + PyArray_Descr *dtypeIn = cluster_dt(); + PyArray_Descr *dtypeOut = cluster_analysis_dt(); + + PyObject *c_obj; + if (!PyArg_ParseTuple(args, "O", &c_obj)) + return NULL; + + // Create two numpy arrays from the passed objects, if possible numpy will + // use the underlying buffer, otherwise it will create a copy, for example + // if data type is different or we pass in a list. The + // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. + PyObject *c_array =PyArray_FromArray((PyArrayObject *)c_obj, dtypeIn, NPY_ARRAY_C_CONTIGUOUS); + + // If parsing of a or b fails we throw an exception in Python + if (c_array == NULL ) { + PyErr_SetString( + PyExc_TypeError, + "Could not convert one of the arguments to a numpy array."); + return NULL; + } + + + const int ndim = PyArray_NDIM((PyArrayObject *)c_array); + + npy_intp *dims = PyArray_SHAPE((PyArrayObject *)c_array); + + Py_ssize_t size=dims[0]; + //printf("%d size %d %d\n",ndim,size,sizeof(ClusterAnalysis)); + //dims[0]=size; + + //Cluster *clusters = reinterpret_cast( PyArray_DATA(reinterpret_cast(c_array))); + + Cluster *clusters = (Cluster *)(PyArray_DATA((PyArrayObject *)(c_array))); + + + + + + // PyObject *PyArray_SimpleNewFromDescr(int nd, npy_int const *dims, PyArray_Descr *descr) + PyObject *clustersA = PyArray_SimpleNewFromDescr(ndim, dims, dtypeOut); + //PyArray_FILLWBYTE((PyArrayObject *)clustersA, 0); //zero initialization can be removed later + npy_intp *strides=PyArray_STRIDES(((PyArrayObject *)(clustersA))); + // printf("strides %d %d\n", strides[0],sizeof(ClusterAnalysis)); + + //Get a pointer to the array memory + ClusterAnalysis* buf = PyArray_DATA((PyArrayObject *)clustersA); + + // Call the standalone C code to read clusters from file + // Here goes the looping, removing frame numbers etc. + int nc = analyze_clusters(size,clusters,buf); + //printf("%d %d\n",nc,size); + + if(nc != size){ + + PyErr_SetString( + PyExc_TypeError, + "Parsed wrong size array!"); + + } + + return clustersA; +} + + + + + // List all methods in our ClusterFileReader class static PyMethodDef ClusterFileReader_methods[] = { {"read", (PyCFunction) ClusterFileReader_read, METH_VARARGS, "Read clusters" }, + {"clusterize", (PyCFunction) ClusterFileReader_clusterize, METH_VARARGS, + "Analyze clusters" + }, {NULL} /* Sentinel */ }; @@ -155,12 +252,14 @@ static PyTypeObject ClusterFileReaderType = { static char module_docstring[] = "C functions to read cluster files"; -// This is the module itself -static PyMethodDef module_methods[] = { - {NULL, NULL, 0, NULL}}; +/* // This is the module itself */ +/* static PyMethodDef module_methods[] = { */ +/* {"read", (PyCFunction)read, METH_VARARGS, add_doc}, {"decode", (PyCFunction)decode, METH_VARARGS, decode_doc}, */ +/* {NULL, NULL, 0, NULL}}; */ + static struct PyModuleDef creader_def = {PyModuleDef_HEAD_INIT, "creader", - module_docstring, -1, module_methods}; + module_docstring, -1, ClusterFileReader_methods}; PyMODINIT_FUNC PyInit_creader(void) { PyObject *m; diff --git a/src/data_types.h b/src/data_types.h index b7e8abd..7a47c1b 100644 --- a/src/data_types.h +++ b/src/data_types.h @@ -9,4 +9,34 @@ typedef struct { int32_t data[9]; } Cluster ; +typedef enum { + cBottomLeft=0, + cBottomRight=1, + cTopLeft=2, + cTopRight=3 +} corner; + +typedef enum { + pBottomLeft=0, + pBottom=1, + pBottomRight=2, + pLeft=3, + pCenter=4, + pRight=5, + pTopLeft=6, + pTop=7, + qTopRight=8 +} pixel; + + +typedef struct { + int32_t tot2; + int32_t tot3; + int8_t c; +} ClusterAnalysis ; + + + + + #endif diff --git a/test.py b/test.py index 4d5c44b..a6cbd4e 100644 --- a/test.py +++ b/test.py @@ -2,13 +2,15 @@ from creader import ClusterFileReader import numpy as np maxph=100000000 -maxph=100 +#maxph=100 from pathlib import Path -fpath = Path('/Users/erik/data/clusters/beam_En700eV_-40deg_300V_10us_d0_f0_100.clust') -# r = ClusterFileReader("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust") +fpath = Path("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust") +# r = ClusterFileReader() r = ClusterFileReader(fpath.as_posix()) a=r.read(maxph) -# print(a[::100000]) +b=r.clusterize(a) +#v=int(maxph/100) +#print(a[::v])