| 99 | |
| 100 | template<typename T> |
| 101 | af_array ccHelper(const Array<T>& img, const Array<uint>& seedx, |
| 102 | const Array<uint>& seedy, const unsigned radius, |
| 103 | const unsigned mult, const unsigned iterations, |
| 104 | const double segmentedValue) { |
| 105 | using CT = |
| 106 | typename conditional<is_same<T, double>::value, double, float>::type; |
| 107 | constexpr CT epsilon = 1.0e-6; |
| 108 | |
| 109 | auto calcVar = [](CT s2, CT s1, CT n) -> CT { |
| 110 | CT retVal = CT(0); |
| 111 | if (n > 1) { retVal = (s2 - (s1 * s1 / n)) / (n - CT(1)); } |
| 112 | return retVal; |
| 113 | }; |
| 114 | |
| 115 | const dim4& inDims = img.dims(); |
| 116 | const dim4& seedDims = seedx.dims(); |
| 117 | const size_t numSeeds = seedx.elements(); |
| 118 | const unsigned nhoodLen = 2 * radius + 1; |
| 119 | const unsigned nhoodSize = nhoodLen * nhoodLen; |
| 120 | |
| 121 | auto labelSegmented = [segmentedValue, inDims](const Array<CT>& segmented) { |
| 122 | Array<CT> newVals = createValueArray(inDims, CT(segmentedValue)); |
| 123 | Array<CT> result = arithOp<CT, af_mul_t>(newVals, segmented, inDims); |
| 124 | // cast final result to input type |
| 125 | return cast<T, CT>(result); |
| 126 | }; |
| 127 | |
| 128 | Array<uint> radiip = createValueArray<uint>(seedDims, radius + 1); |
| 129 | Array<uint> radii = createValueArray<uint>(seedDims, radius); |
| 130 | Array<uint> _x = arithOp<uint, af_sub_t>(seedx, radiip, seedDims); |
| 131 | Array<uint> x_ = arithOp<uint, af_add_t>(seedx, radii, seedDims); |
| 132 | Array<uint> _y = arithOp<uint, af_sub_t>(seedy, radiip, seedDims); |
| 133 | Array<uint> y_ = arithOp<uint, af_add_t>(seedy, radii, seedDims); |
| 134 | Array<CT> in = convRange<CT, T>(img, CT(1), CT(2)); |
| 135 | Array<CT> in_2 = arithOp<CT, af_mul_t>(in, in, inDims); |
| 136 | Array<CT> I1 = integralImage<CT>(in); |
| 137 | Array<CT> I2 = integralImage<CT>(in_2); |
| 138 | Array<CT> S1 = sum(I1, _x, x_, _y, y_); |
| 139 | Array<CT> S2 = sum(I2, _x, x_, _y, y_); |
| 140 | CT totSum = getScalar<CT>(reduce_all<af_add_t, CT, CT>(S1)); |
| 141 | CT totSumSq = getScalar<CT>(reduce_all<af_add_t, CT, CT>(S2)); |
| 142 | CT totalNum = numSeeds * nhoodSize; |
| 143 | CT s1mean = totSum / totalNum; |
| 144 | CT s1var = calcVar(totSumSq, totSum, totalNum); |
| 145 | CT s1stddev = sqrt(s1var); |
| 146 | CT lower = s1mean - mult * s1stddev; |
| 147 | CT upper = s1mean + mult * s1stddev; |
| 148 | |
| 149 | Array<CT> seedIntensities = pointList(in, seedx, seedy); |
| 150 | CT maxSeedIntensity = |
| 151 | getScalar<CT>(reduce_all<af_max_t, CT, CT>(seedIntensities)); |
| 152 | CT minSeedIntensity = |
| 153 | getScalar<CT>(reduce_all<af_min_t, CT, CT>(seedIntensities)); |
| 154 | |
| 155 | if (lower > minSeedIntensity) { lower = minSeedIntensity; } |
| 156 | if (upper < maxSeedIntensity) { upper = maxSeedIntensity; } |
| 157 | |
| 158 | Array<CT> segmented = floodFill(in, seedx, seedy, CT(1), lower, upper); |
no test coverage detected