Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	src/pull/ccPull.c	Lines:	0	123	0.0 %
Date:	2017-05-26	Branches:	0	80	0.0 %


/*
  Teem: Tools to process and visualize scientific data and images             .
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library 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
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/

#include "pull.h"
#include "privatePull.h"

/*
** HEY: this messes with the points' idtag, and pctx->idtagNext,
** even though it really shouldn't have to
*/
int
pullCCFind(pullContext *pctx) {
  static const char me[]="pullCCFind";
  airArray *mop, *eqvArr;
  unsigned int passIdx, binIdx, pointIdx, neighIdx, eqvNum,
    pointNum, *idmap;
  pullBin *bin;
  pullPoint *point, *her;

  if (!pctx) {
    biffAddf(PULL, "%s: got NULL pointer", me);
    return 1;
  }
  if (_pullIterate(pctx, pullProcessModeNeighLearn)) {
    biffAddf(PULL, "%s: trouble with %s for CC", me,
             airEnumStr(pullProcessMode, pullProcessModeNeighLearn));
    return 1;
  }

  mop = airMopNew();
  pointNum = pullPointNumber(pctx);
  eqvArr = airArrayNew(NULL, NULL, 2*sizeof(unsigned int), pointNum);
  airMopAdd(mop, eqvArr, (airMopper)airArrayNuke, airMopAlways);
  idmap = AIR_CAST(unsigned int *, calloc(pointNum, sizeof(unsigned int)));
  airMopAdd(mop, idmap, airFree, airMopAlways);

  /* to be safe, renumber all points, so that we know that the
     idtags are contiguous, starting at 0. HEY: this should handled
     by having a map from real point->idtag to a point number assigned
     just for the sake of doing CCs */
  pctx->idtagNext = 0;
  for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
    bin = pctx->bin + binIdx;
    for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
      point = bin->point[pointIdx];
      point->idtag = pctx->idtagNext++;
    }
  }

  /* same stupidity copied from limn/polymod.c:limnPolyDataCCFind */
  eqvNum = 0;
  for (passIdx=0; passIdx<2; passIdx++) {
    if (passIdx) {
      airArrayLenPreSet(eqvArr, eqvNum);
    }
    for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
      bin = pctx->bin + binIdx;
      for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
        point = bin->point[pointIdx];
        for (neighIdx=0; neighIdx<point->neighPointNum; neighIdx++) {
          if (0 == passIdx) {
            ++eqvNum;
          } else {
            her = point->neighPoint[neighIdx];
            airEqvAdd(eqvArr, point->idtag, her->idtag);
          }
        }
      }
    }
  }

  /* do the CC analysis */
  pctx->CCNum = airEqvMap(eqvArr, idmap, pointNum);

  /* assign idcc's */
  for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
    bin = pctx->bin + binIdx;
    for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
      point = bin->point[pointIdx];
      point->idCC = idmap[point->idtag];
    }
  }

  airMopOkay(mop);
  return 0;
}

/*
** measure something about connected componts already found
**
** measrInfo can be 0 to say "measure # particles in CC", or
** it can be a scalar pullInfo
**
** rho == 0: only size, rho == 1: only measrInfo
*/
int
pullCCMeasure(pullContext *pctx, Nrrd *nmeasr, int measrInfo, double rho) {
  static const char me[]="pullCCMeasure";
  airArray *mop;
  unsigned int binIdx, pointIdx, ii;
  double *meas, *size;
  pullBin *bin;
  pullPoint *point;

  if (!( pctx && nmeasr )) {
    biffAddf(PULL, "%s: got NULL pointer", me);
    return 1;
  }
  if (!pctx->CCNum) {
    biffAddf(PULL, "%s: CCNum == 0: haven't yet learned CCs?", me);
    return 1;
  }
  if (measrInfo) {
    if (airEnumValCheck(pullInfo, measrInfo)) {
      biffAddf(PULL, "%s: measrInfo %d not a valid %s", me,
               measrInfo, pullInfo->name);
      return 1;
    }
    if (1 != pullInfoLen(measrInfo)) {
      biffAddf(PULL, "%s: measrInfo %s (%d) isn't a scalar (len %u)", me,
               airEnumStr(pullInfo, measrInfo), measrInfo,
               pullInfoLen(measrInfo));
      return 1;
    }
    if (!pctx->ispec[measrInfo]) {
      biffAddf(PULL, "%s: no ispec set for measrInfo %s (%d)", me,
               airEnumStr(pullInfo, measrInfo), measrInfo);
      return 1;
    }
  } /* else measrInfo is zero, they want to know # points */
  /* in any case nmeasr is allocated for doubles */
  if (nrrdMaybeAlloc_va(nmeasr, nrrdTypeDouble, 1,
                        AIR_CAST(size_t, pctx->CCNum))) {
    biffMovef(PULL, NRRD, "%s: couldn't alloc nmeasr", me);
    return 1;
  }
  meas = AIR_CAST(double *, nmeasr->data);

  mop = airMopNew();
  /* HEY: don't actually need to allocate and set size[],
     if measrInfo == 0 */
  if (!(size = AIR_CAST(double *,
                        calloc(pctx->CCNum, sizeof(double))))) {
    biffAddf(PULL, "%s: couldn't alloc size", me);
    airMopError(mop); return 1;
  }
  airMopAdd(mop, size, airFree, airMopAlways);

  /* finally, do measurement */
  for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
    bin = pctx->bin + binIdx;
    for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
      point = bin->point[pointIdx];
      size[point->idCC]++;
      meas[point->idCC] += (measrInfo
                            ? pullPointScalar(pctx, point, measrInfo,
                                              NULL, NULL)
                            : 1);
    }
  }
  if (measrInfo) {
    for (ii=0; ii<pctx->CCNum; ii++) {
      meas[ii] /= size[ii];
      meas[ii] = AIR_LERP(rho, size[ii], meas[ii]);
    }
  }

  airMopOkay(mop);
  return 0;
}

typedef struct {
  unsigned int i;
  double d;
} ccpair;

/* we intend to sort in *descending* order */
static int
ccpairCompare(const void *_a, const void *_b) {
  const ccpair *a, *b;

  a = AIR_CAST(const ccpair *, _a);
  b = AIR_CAST(const ccpair *, _b);
  return (a->d < b->d
          ? +1
          : (a->d > b->d
             ? -1
             : 0));
}

/*
******** pullCCSort
**
** sort CCs in pull context
**
** measrInfo == 0: sort by size, else,
** sort by blend of size and measrInfo
** rho == 0: only size, rho == 1: only measrInfo
*/
int
pullCCSort(pullContext *pctx, int measrInfo, double rho) {
  static const char me[]="pullCCSort";
  ccpair *pair;
  Nrrd *nmeasr;
  airArray *mop;
  unsigned int ii, *revm, binIdx, pointIdx;
  double *measr;
  pullBin *bin;
  pullPoint *point;
  int E;

  if (!pctx) {
    biffAddf(PULL, "%s: got NULL pointer", me);
    return 1;
  }
  if (!pctx->CCNum) {
    biffAddf(PULL, "%s: haven't yet learned CCs?", me);
    return 1;
  }

#define CALLOC(T) (AIR_CAST(T *, calloc(pctx->CCNum, sizeof(T))))
  mop = airMopNew();
  if (!(nmeasr = nrrdNew())
      || airMopAdd(mop, nmeasr, (airMopper)nrrdNuke, airMopAlways)
      || !(pair = CALLOC(ccpair))
      || airMopAdd(mop, pair, airFree, airMopAlways)
      || !(revm = CALLOC(unsigned int))
      || airMopAdd(mop, revm, airFree, airMopAlways)) {
    biffAddf(PULL, "%s: couldn't allocate everything", me);
    airMopError(mop); return 1;
  }
#undef CALLOC
  if (!measrInfo) {
    /* sorting by size */
    E = pullCCMeasure(pctx, nmeasr, 0, 0.0);
  } else {
    /* sorting by some blend of size and pullInfo */
    E = pullCCMeasure(pctx, nmeasr, measrInfo, rho);
  }
  if (E) {
    biffAddf(PULL, "%s: problem measuring CCs", me);
    airMopError(mop); return 1;
  }

  measr = AIR_CAST(double *, nmeasr->data);
  for (ii=0; ii<pctx->CCNum; ii++) {
    pair[ii].i = ii;
    pair[ii].d = measr[ii];
  }
  qsort(pair, pctx->CCNum, sizeof(ccpair), ccpairCompare);
  for (ii=0; ii<pctx->CCNum; ii++) {
    revm[pair[ii].i] = ii;
  }

  for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
    bin = pctx->bin + binIdx;
    for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
      point = bin->point[pointIdx];
      point->idCC = revm[point->idCC];
    }
  }

  airMopOkay(mop);
  return 0;
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/*
2			Teem: Tools to process and visualize scientific data and images .
3			Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4			Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5			Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7			This library is free software; you can redistribute it and/or
8			modify it under the terms of the GNU Lesser General Public License
9			(LGPL) as published by the Free Software Foundation; either
10			version 2.1 of the License, or (at your option) any later version.
11			The terms of redistributing and/or modifying this software also
12			include exceptions to the LGPL that facilitate static linking.
13
14			This library is distributed in the hope that it will be useful,
15			but WITHOUT ANY WARRANTY; without even the implied warranty of
16			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17			Lesser General Public License for more details.
18
19			You should have received a copy of the GNU Lesser General Public License
20			along with this library; if not, write to Free Software Foundation, Inc.,
21			51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22			*/
23
24			#include "pull.h"
25			#include "privatePull.h"
26
27			/*
28			** HEY: this messes with the points' idtag, and pctx->idtagNext,
29			** even though it really shouldn't have to
30			*/
31			int
32			pullCCFind(pullContext *pctx) {
33			static const char me[]="pullCCFind";
34			airArray mop, eqvArr;
35			unsigned int passIdx, binIdx, pointIdx, neighIdx, eqvNum,
36			pointNum, *idmap;
37			pullBin *bin;
38			pullPoint point, her;
39
40			if (!pctx) {
41			biffAddf(PULL, "%s: got NULL pointer", me);
42			return 1;
43			}
44			if (_pullIterate(pctx, pullProcessModeNeighLearn)) {
45			biffAddf(PULL, "%s: trouble with %s for CC", me,
46			airEnumStr(pullProcessMode, pullProcessModeNeighLearn));
47			return 1;
48			}
49
50			mop = airMopNew();
51			pointNum = pullPointNumber(pctx);
52			eqvArr = airArrayNew(NULL, NULL, 2*sizeof(unsigned int), pointNum);
53			airMopAdd(mop, eqvArr, (airMopper)airArrayNuke, airMopAlways);
54			idmap = AIR_CAST(unsigned int *, calloc(pointNum, sizeof(unsigned int)));
55			airMopAdd(mop, idmap, airFree, airMopAlways);
56
57			/* to be safe, renumber all points, so that we know that the
58			idtags are contiguous, starting at 0. HEY: this should handled
59			by having a map from real point->idtag to a point number assigned
60			just for the sake of doing CCs */
61			pctx->idtagNext = 0;
62			for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
63			bin = pctx->bin + binIdx;
64			for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
65			point = bin->point[pointIdx];
66			point->idtag = pctx->idtagNext++;
67			}
68			}
69
70			/* same stupidity copied from limn/polymod.c:limnPolyDataCCFind */
71			eqvNum = 0;
72			for (passIdx=0; passIdx<2; passIdx++) {
73			if (passIdx) {
74			airArrayLenPreSet(eqvArr, eqvNum);
75			}
76			for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
77			bin = pctx->bin + binIdx;
78			for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
79			point = bin->point[pointIdx];
80			for (neighIdx=0; neighIdx<point->neighPointNum; neighIdx++) {
81			if (0 == passIdx) {
82			++eqvNum;
83			} else {
84			her = point->neighPoint[neighIdx];
85			airEqvAdd(eqvArr, point->idtag, her->idtag);
86			}
87			}
88			}
89			}
90			}
91
92			/* do the CC analysis */
93			pctx->CCNum = airEqvMap(eqvArr, idmap, pointNum);
94
95			/* assign idcc's */
96			for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
97			bin = pctx->bin + binIdx;
98			for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
99			point = bin->point[pointIdx];
100			point->idCC = idmap[point->idtag];
101			}
102			}
103
104			airMopOkay(mop);
105			return 0;
106			}
107
108			/*
109			** measure something about connected componts already found
110			**
111			** measrInfo can be 0 to say "measure # particles in CC", or
112			** it can be a scalar pullInfo
113			**
114			** rho == 0: only size, rho == 1: only measrInfo
115			*/
116			int
117			pullCCMeasure(pullContext pctx, Nrrd nmeasr, int measrInfo, double rho) {
118			static const char me[]="pullCCMeasure";
119			airArray *mop;
120			unsigned int binIdx, pointIdx, ii;
121			double meas, size;
122			pullBin *bin;
123			pullPoint *point;
124
125			if (!( pctx && nmeasr )) {
126			biffAddf(PULL, "%s: got NULL pointer", me);
127			return 1;
128			}
129			if (!pctx->CCNum) {
130			biffAddf(PULL, "%s: CCNum == 0: haven't yet learned CCs?", me);
131			return 1;
132			}
133			if (measrInfo) {
134			if (airEnumValCheck(pullInfo, measrInfo)) {
135			biffAddf(PULL, "%s: measrInfo %d not a valid %s", me,
136			measrInfo, pullInfo->name);
137			return 1;
138			}
139			if (1 != pullInfoLen(measrInfo)) {
140			biffAddf(PULL, "%s: measrInfo %s (%d) isn't a scalar (len %u)", me,
141			airEnumStr(pullInfo, measrInfo), measrInfo,
142			pullInfoLen(measrInfo));
143			return 1;
144			}
145			if (!pctx->ispec[measrInfo]) {
146			biffAddf(PULL, "%s: no ispec set for measrInfo %s (%d)", me,
147			airEnumStr(pullInfo, measrInfo), measrInfo);
148			return 1;
149			}
150			} /* else measrInfo is zero, they want to know # points */
151			/* in any case nmeasr is allocated for doubles */
152			if (nrrdMaybeAlloc_va(nmeasr, nrrdTypeDouble, 1,
153			AIR_CAST(size_t, pctx->CCNum))) {
154			biffMovef(PULL, NRRD, "%s: couldn't alloc nmeasr", me);
155			return 1;
156			}
157			meas = AIR_CAST(double *, nmeasr->data);
158
159			mop = airMopNew();
160			/* HEY: don't actually need to allocate and set size[],
161			if measrInfo == 0 */
162			if (!(size = AIR_CAST(double *,
163			calloc(pctx->CCNum, sizeof(double))))) {
164			biffAddf(PULL, "%s: couldn't alloc size", me);
165			airMopError(mop); return 1;
166			}
167			airMopAdd(mop, size, airFree, airMopAlways);
168
169			/* finally, do measurement */
170			for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
171			bin = pctx->bin + binIdx;
172			for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
173			point = bin->point[pointIdx];
174			size[point->idCC]++;
175			meas[point->idCC] += (measrInfo
176			? pullPointScalar(pctx, point, measrInfo,
177			NULL, NULL)
178			: 1);
179			}
180			}
181			if (measrInfo) {
182			for (ii=0; ii<pctx->CCNum; ii++) {
183			meas[ii] /= size[ii];
184			meas[ii] = AIR_LERP(rho, size[ii], meas[ii]);
185			}
186			}
187
188			airMopOkay(mop);
189			return 0;
190			}
191
192			typedef struct {
193			unsigned int i;
194			double d;
195			} ccpair;
196
197			/* we intend to sort in descending order */
198			static int
199			ccpairCompare(const void _a, const void _b) {
200			const ccpair a, b;
201
202			a = AIR_CAST(const ccpair *, _a);
203			b = AIR_CAST(const ccpair *, _b);
204			return (a->d < b->d
205			? +1
206			: (a->d > b->d
207			? -1
208			: 0));
209			}
210
211			/*
212			******** pullCCSort
213			**
214			** sort CCs in pull context
215			**
216			** measrInfo == 0: sort by size, else,
217			** sort by blend of size and measrInfo
218			** rho == 0: only size, rho == 1: only measrInfo
219			*/
220			int
221			pullCCSort(pullContext *pctx, int measrInfo, double rho) {
222			static const char me[]="pullCCSort";
223			ccpair *pair;
224			Nrrd *nmeasr;
225			airArray *mop;
226			unsigned int ii, *revm, binIdx, pointIdx;
227			double *measr;
228			pullBin *bin;
229			pullPoint *point;
230			int E;
231
232			if (!pctx) {
233			biffAddf(PULL, "%s: got NULL pointer", me);
234			return 1;
235			}
236			if (!pctx->CCNum) {
237			biffAddf(PULL, "%s: haven't yet learned CCs?", me);
238			return 1;
239			}
240
241			#define CALLOC(T) (AIR_CAST(T *, calloc(pctx->CCNum, sizeof(T))))
242			mop = airMopNew();
243			if (!(nmeasr = nrrdNew())
244			\|\| airMopAdd(mop, nmeasr, (airMopper)nrrdNuke, airMopAlways)
245			\|\| !(pair = CALLOC(ccpair))
246			\|\| airMopAdd(mop, pair, airFree, airMopAlways)
247			\|\| !(revm = CALLOC(unsigned int))
248			\|\| airMopAdd(mop, revm, airFree, airMopAlways)) {
249			biffAddf(PULL, "%s: couldn't allocate everything", me);
250			airMopError(mop); return 1;
251			}
252			#undef CALLOC
253			if (!measrInfo) {
254			/* sorting by size */
255			E = pullCCMeasure(pctx, nmeasr, 0, 0.0);
256			} else {
257			/* sorting by some blend of size and pullInfo */
258			E = pullCCMeasure(pctx, nmeasr, measrInfo, rho);
259			}
260			if (E) {
261			biffAddf(PULL, "%s: problem measuring CCs", me);
262			airMopError(mop); return 1;
263			}
264
265			measr = AIR_CAST(double *, nmeasr->data);
266			for (ii=0; ii<pctx->CCNum; ii++) {
267			pair[ii].i = ii;
268			pair[ii].d = measr[ii];
269			}
270			qsort(pair, pctx->CCNum, sizeof(ccpair), ccpairCompare);
271			for (ii=0; ii<pctx->CCNum; ii++) {
272			revm[pair[ii].i] = ii;
273			}
274
275			for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
276			bin = pctx->bin + binIdx;
277			for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
278			point = bin->point[pointIdx];
279			point->idCC = revm[point->idCC];
280			}
281			}
282
283			airMopOkay(mop);
284			return 0;
285			}