1.1 --- a/src/mask.h Thu Jun 11 15:20:52 2009 +0000
1.2 +++ b/src/mask.h Sun Jun 28 01:26:46 2009 -0700
1.3 @@ -347,11 +347,24 @@
1.4 // mem usage after: CoarseMask: 2/8 * uBB * MaskType
1.5 // !CoarseMask: 2 * uBB * MaskType
1.6
1.7 +ImageExportInfo nftmaskinfo1("enblend_nft_input.tif");
1.8 +exportImage(srcImageRange(*nftInputImage), nftmaskinfo1);
1.9 +
1.10 nearestFeatureTransform(wraparound,
1.11 srcImageRange(*nftInputImage),
1.12 destIter(nftOutputImage->upperLeft() + nftOutputOffset),
1.13 NumericTraits<MaskPixelType>::one());
1.14
1.15 +ImageExportInfo nftmaskinfo2("enblend_nft_output.tif");
1.16 +exportImage(srcImageRange(*nftOutputImage), nftmaskinfo2);
1.17 +
1.18 + nearestFeatureTransformExact(srcImageRange(*nftInputImage),
1.19 + destIter(nftOutputImage->upperLeft() + nftOutputOffset),
1.20 + NumericTraits<MaskPixelType>::one());
1.21 +
1.22 +ImageExportInfo nftmaskinfo3("enblend_nfte_output.tif");
1.23 +exportImage(srcImageRange(*nftOutputImage), nftmaskinfo3);
1.24 +
1.25 //ImageExportInfo nftmaskinfo("enblend_nft.tif");
1.26 //exportImage(srcImageRange(*nftOutputImage), nftmaskinfo);
1.27
2.1 --- a/src/nearest.h Thu Jun 11 15:20:52 2009 +0000
2.2 +++ b/src/nearest.h Sun Jun 28 01:26:46 2009 -0700
2.3 @@ -167,6 +167,318 @@
2.4
2.5 };
2.6
2.7 +template<class FeatureType> class Chain {
2.8 +public:
2.9 + struct Node {
2.10 + bool active;
2.11 + FeatureType f;
2.12 + int x;
2.13 + int y;
2.14 + Node* left;
2.15 + Node* right;
2.16 + };
2.17 +
2.18 + Chain(int w) {
2.19 + nodeArray = new Node[w];
2.20 + activeNodes = 0;
2.21 + leftmostNode = NULL;
2.22 + rightmostNode = NULL;
2.23 + for (int i = 0; i < w; ++i) {
2.24 + nodeArray[i].active = false;
2.25 + nodeArray[i].x = i;
2.26 + nodeArray[i].left = NULL;
2.27 + nodeArray[i].right = NULL;
2.28 + }
2.29 + }
2.30 +
2.31 + ~Chain() {
2.32 + delete [] nodeArray;
2.33 + }
2.34 +
2.35 + void insertAfter(Node* n, int x, int y, FeatureType f) {
2.36 + if (nodeArray[x].active) {
2.37 + nodeArray[x].y = y;
2.38 + nodeArray[x].f = f;
2.39 + }
2.40 + else if (n == NULL) {
2.41 + // Insert at beginning
2.42 + nodeArray[x].active = true;
2.43 + nodeArray[x].f = f;
2.44 + nodeArray[x].y = y;
2.45 + nodeArray[x].left = NULL;
2.46 + nodeArray[x].right = leftmostNode;
2.47 + if (leftmostNode != NULL) {
2.48 + leftmostNode->left = &nodeArray[x];
2.49 + }
2.50 + leftmostNode = &nodeArray[x];
2.51 + if (rightmostNode == NULL) {
2.52 + rightmostNode = &nodeArray[x];
2.53 + }
2.54 + ++activeNodes;
2.55 + } else {
2.56 + // Insert after given node.
2.57 + nodeArray[x].active = true;
2.58 + nodeArray[x].f = f;
2.59 + nodeArray[x].y = y;
2.60 + nodeArray[x].left = n;
2.61 + nodeArray[x].right = n->right;
2.62 + if (n->right == NULL) {
2.63 + rightmostNode = &nodeArray[x];
2.64 + } else {
2.65 + nodeArray[x].right->left = &nodeArray[x];
2.66 + }
2.67 + n->right = &nodeArray[x];
2.68 + ++activeNodes;
2.69 + }
2.70 + }
2.71 +
2.72 + void erase(Node* n) {
2.73 + n->active = false;
2.74 + if (n->left == NULL) {
2.75 + leftmostNode = n->right;
2.76 + } else {
2.77 + n->left->right = n->right;
2.78 + }
2.79 + if (n->right == NULL) {
2.80 + rightmostNode = n->left;
2.81 + } else {
2.82 + n->right->left = n->left;
2.83 + }
2.84 + n->left = NULL;
2.85 + n->right = NULL;
2.86 + --activeNodes;
2.87 + }
2.88 +
2.89 + // Compute the y-coordinate of the intersection of the
2.90 + // perpendicular bisector of pq with the vertical line at x.
2.91 + pair<double, bool> pbY(Node* p, double x) {
2.92 + Node* q = p->right;
2.93 + assert(q != NULL);
2.94 + const double px = p->x;
2.95 + const double py = p->y;
2.96 + const double qx = q->x;
2.97 + const double qy = q->y;
2.98 + if (q->y == p->y) return std::make_pair<double, bool>(0.0, false);
2.99 + const double midX = (qx + px) / 2.0;
2.100 + const double midY = (qy + py) / 2.0;
2.101 + const double m = (qy - py) / (qx - px);
2.102 + const double mPB = -1.0 / m;
2.103 + // y = mPB * (x - midX) + midY
2.104 + return std::make_pair<double, bool>((mPB * (x - midX) + midY), true);
2.105 + }
2.106 +
2.107 + // Compute the x-coordinate of the intersection of the
2.108 + // perpendicular bisector of pq with the horizontal line at y.
2.109 + double pbX(Node* p, double y) {
2.110 + Node* q = p->right;
2.111 + assert(q != NULL);
2.112 + const double px = p->x;
2.113 + const double py = p->y;
2.114 + const double qx = q->x;
2.115 + const double qy = q->y;
2.116 + const double midX = (qx + px) / 2.0;
2.117 + if (q->y == p->y) return midX;
2.118 + const double midY = (qy + py) / 2.0;
2.119 + const double m = (qy - py) / (qx - px);
2.120 + const double mPB = -1.0 / m;
2.121 + // y = mPB * (x - midX) + midY
2.122 + // y - midY = mPB * (x - midX)
2.123 + // (y - midY) / mPB = x - midX
2.124 + // x = (y - midY) / mPB + midX
2.125 + return ((y - midY) / mPB + midX);
2.126 + }
2.127 +
2.128 + // Compute the y-coordinate of the intersection of the
2.129 + // perpendicular bisectors of pq and qr.
2.130 + pair<double, bool> intersectionY(Node* p) {
2.131 + Node* q = p->right;
2.132 + assert(q != NULL);
2.133 + Node* r = q->right;
2.134 + assert(r != NULL);
2.135 + const int px = p->x;
2.136 + const int py = p->y;
2.137 + const int qx = q->x;
2.138 + const int qy = q->y;
2.139 + const int rx = r->x;
2.140 + const int ry = r->y;
2.141 +
2.142 + const double midPQX = (qx + px) / 2.0;
2.143 + if (qy == py) return pbY(q, midPQX);
2.144 +
2.145 + const double midQRX = (rx + qx) / 2.0;
2.146 + if (qy == ry) return pbY(p, midQRX);
2.147 +
2.148 + const double midPQY = (qy + py) / 2.0;
2.149 + const double slopePQ = static_cast<double>(qy - py) / static_cast<double>(qx - px);
2.150 + const double pbSlopePQ = -1.0 / slopePQ;
2.151 +
2.152 + const double midQRY = (ry + qy) / 2.0;
2.153 + const double slopeQR = static_cast<double>(ry - qy) / static_cast<double>(rx - qx);
2.154 + const double pbSlopeQR = -1.0 / slopeQR;
2.155 +
2.156 + if (slopePQ == slopeQR) return std::make_pair<double, bool>(0.0, false);
2.157 +
2.158 + // pbPQ: y = pbSlopePQ * (x - midPQX) + midPQY;
2.159 + // pbQR: y = pbSlopeQR * (x - midQRX) + midQRY;
2.160 + // y = (pbSlopeQR * x) - (pbSlopeQR * midQRX) + midQRY
2.161 + // pbSlopeQR * x = y + (pbSlopeQR * midQRX) - midQRY
2.162 + // x = (y + (pbSlopeQR * midQRX) - midQRY) / pbSlopeQR
2.163 + // substituting into pbPQ:
2.164 + // y = pbSlopePQ * (y + (pbSlopeQR * midQRX) - midQRY) / pbSlopeQR - pbSlopePQ*midPQX + midPQY;
2.165 + // = (pbSlopePQ/pbSlopeQR) * (y + pbSlopeQR*midQRX - midQRY) - pbSlopePQ*midPQX + midPQY;
2.166 + // y - (pbSlopePQ/pbSlopeQR)*y = (pbSlopePQ/pbSlopeQR)*(pbSlopeQR*midQRX - midQRY) - pbSlopePQ*midPQX + midPQY
2.167 + // y * (1 - pbSlopePQ/pbSlopeQR) = ...
2.168 + // y = ... / (1 - pbSlopePQ/pbSlopeQR)
2.169 + const double iy = ((pbSlopePQ/pbSlopeQR)*(pbSlopeQR*midQRX - midQRY) - pbSlopePQ*midPQX + midPQY)
2.170 + / (1.0 - pbSlopePQ / pbSlopeQR);
2.171 + return std::make_pair<double, bool>(iy, true);
2.172 + }
2.173 +
2.174 + Node* operator[](int idx) {
2.175 + return &(nodeArray[idx]);
2.176 + }
2.177 +
2.178 + Node* nodeArray;
2.179 + int activeNodes;
2.180 + Node* leftmostNode;
2.181 + Node* rightmostNode;
2.182 +};
2.183 +
2.184 +template <class SrcImageIterator, class SrcAccessor,
2.185 + class DestImageIterator, class DestAccessor>
2.186 +void nearestFeatureTransformExact(
2.187 + SrcImageIterator src_upperleft,
2.188 + SrcImageIterator src_lowerright,
2.189 + SrcAccessor sa,
2.190 + DestImageIterator dest_upperleft,
2.191 + DestAccessor da,
2.192 + typename SrcAccessor::value_type featurelessPixel) {
2.193 +
2.194 + typedef typename EnblendNumericTraits<UInt32>::ImageType DNFImage;
2.195 + typedef typename DNFImage::traverser DnfIterator;
2.196 + typedef typename SrcAccessor::value_type SrcValueType;
2.197 +
2.198 + SrcImageIterator sx, sy, send;
2.199 + DnfIterator dnfAboveX, dnfAboveY;
2.200 + DestImageIterator dx, dy;
2.201 +
2.202 + int w = src_lowerright.x - src_upperleft.x;
2.203 + int h = src_lowerright.y - src_upperleft.y;
2.204 +
2.205 + DNFImage dnfAbove(w, h, NumericTraits<UInt32>::max());
2.206 + Chain<SrcValueType> chain(w);
2.207 +
2.208 + // Top-down pass
2.209 + sy = src_upperleft;
2.210 + send = src_lowerright;
2.211 + dy = dest_upperleft;
2.212 + dnfAboveY = dnfAbove.upperLeft();
2.213 +
2.214 + for (int yIndex = 0; sy.y < send.y; ++yIndex, ++sy.y, ++dy.y, ++dnfAboveY.y) {
2.215 + // Insert new features into chain.
2.216 + sx = sy;
2.217 + typename Chain<SrcValueType>::Node* lastActiveNode = NULL;
2.218 + for (int xIndex = 0; sx.x < send.x; ++xIndex, ++sx.x) {
2.219 + SrcValueType pixel = sa(sx);
2.220 + if (pixel != featurelessPixel) {
2.221 + //cout << "new feature xIndex=" << xIndex << " yIndex=" << yIndex << " f=" << (int)pixel << endl;
2.222 + chain.insertAfter(lastActiveNode, xIndex, yIndex, pixel);
2.223 + }
2.224 + if (chain[xIndex]->active) lastActiveNode = chain[xIndex];
2.225 + }
2.226 +
2.227 + //cout << "after insert: ";
2.228 + //for (lastActiveNode = chain.leftmostNode; lastActiveNode != NULL; lastActiveNode = lastActiveNode->right) {
2.229 + // cout << "[" << lastActiveNode->x << "," << lastActiveNode->y << ": " << (int)lastActiveNode->f << "] ";
2.230 + //}
2.231 + //cout << endl;
2.232 +
2.233 + // Remove leftmost extreme features
2.234 + while (chain.activeNodes > 1) {
2.235 + // If intersection of perpendicular bisector of pq and x=0 is < yIndex, delete p.
2.236 + pair<double, bool> r = chain.pbY(chain.leftmostNode, 0);
2.237 + if (r.second && (r.first < yIndex) && (r.first > chain.leftmostNode->y)) chain.erase(chain.leftmostNode);
2.238 + else break;
2.239 + }
2.240 +
2.241 + //cout << "after leftmost: ";
2.242 + //for (lastActiveNode = chain.leftmostNode; lastActiveNode != NULL; lastActiveNode = lastActiveNode->right) {
2.243 + // cout << "[" << lastActiveNode->x << "," << lastActiveNode->y << ": " << (int)lastActiveNode->f << "] ";
2.244 + //}
2.245 +
2.246 + // Remove rightmost extreme features
2.247 + while (chain.activeNodes > 1) {
2.248 + // If intersection of perpendicular bisector of pq and x=w-1 is <= yIndex, delete q.
2.249 + pair<double, bool> r = chain.pbY(chain.rightmostNode->left, w-1);
2.250 + if (r.second && (r.first <= yIndex) && (r.first > chain.rightmostNode->y)) chain.erase(chain.rightmostNode);
2.251 + else break;
2.252 + }
2.253 +
2.254 + //cout << "after rightmost: ";
2.255 + //for (lastActiveNode = chain.leftmostNode; lastActiveNode != NULL; lastActiveNode = lastActiveNode->right) {
2.256 + // cout << "[" << lastActiveNode->x << "," << lastActiveNode->y << ": " << (int)lastActiveNode->f << "] ";
2.257 + //}
2.258 +
2.259 + // Remove internal features
2.260 + typename Chain<SrcValueType>::Node* p = chain.leftmostNode;
2.261 + while (chain.activeNodes > 2) {
2.262 + typename Chain<SrcValueType>::Node* q = p->right;
2.263 + typename Chain<SrcValueType>::Node* r = q->right;
2.264 + if (r == NULL) break;
2.265 + pair<double, bool> iy = chain.intersectionY(p);
2.266 + if (iy.second && (iy.first < yIndex) && (iy.first > q->y)) {
2.267 + // remove q and backtrack
2.268 + chain.erase(q);
2.269 + if (p->left != NULL) p = p->left;
2.270 + } else {
2.271 + // advance
2.272 + p = q;
2.273 + }
2.274 + }
2.275 +
2.276 + //cout << "after internal: ";
2.277 + //for (lastActiveNode = chain.leftmostNode; lastActiveNode != NULL; lastActiveNode = lastActiveNode->right) {
2.278 + // cout << "[" << lastActiveNode->x << "," << lastActiveNode->y << ": " << (int)lastActiveNode->f << "] ";
2.279 + //}
2.280 + //cout << endl;
2.281 +
2.282 + // assign
2.283 + dx = dy;
2.284 + dnfAboveX = dnfAboveY;
2.285 + lastActiveNode = chain.leftmostNode;
2.286 + int xIndex = 0;
2.287 + while (lastActiveNode != NULL) {
2.288 + int xLimit = (lastActiveNode->right) ? static_cast<int>(floor(chain.pbX(lastActiveNode, yIndex))) : w;
2.289 + //cout << "[" << lastActiveNode->x << "," << lastActiveNode->y << ": " << (int)lastActiveNode->f << "] xLimit=" << xLimit << endl;
2.290 + for (; xIndex < w && xIndex <= xLimit; ++xIndex) {
2.291 + dx(xIndex, 0) = lastActiveNode->f;
2.292 + dnfAboveX(xIndex, 0) = (xIndex - lastActiveNode->x) * (xIndex - lastActiveNode->x) + (yIndex - lastActiveNode->y) * (yIndex - lastActiveNode->y);
2.293 + }
2.294 + lastActiveNode = lastActiveNode->right;
2.295 + }
2.296 +
2.297 + } /* end of top-down pass */
2.298 +/*
2.299 + sy = src_lowerright;
2.300 + send = src_upperleft;
2.301 + dy = dest_upperleft + Diff2D(w, h);
2.302 + dnfAboveY = dnfAbove.lowerRight();
2.303 +
2.304 + for (int yIndex = h; sy.y > send.y;) {
2.305 + --yIndex;
2.306 + --sy.y;
2.307 + --dy.y;
2.308 + --dnfAboveY.y;
2.309 +
2.310 + sx = sy;
2.311 + typename Chain<SrcValueType>::Node* lastActiveNode = NULL;
2.312 + for (int xIndex = w; sx.x > send.x;) {
2.313 + --xIndex;
2.314 + --sx.x;
2.315 +*/
2.316 +
2.317 +}
2.318 +
2.319 /** Compute the nearest feature transform.
2.320 * A non-zero pixel in the src image is considered a feature.
2.321 * Each pixel in the dest image is given the value of the nearest feature
2.322 @@ -550,6 +862,19 @@
2.323
2.324 };
2.325
2.326 +template <class SrcImageIterator, class SrcAccessor,
2.327 + class DestImageIterator, class DestAccessor>
2.328 +inline void nearestFeatureTransformExact(
2.329 + triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2.330 + pair<DestImageIterator, DestAccessor> dest,
2.331 + typename SrcAccessor::value_type featurelessPixel) {
2.332 +
2.333 + nearestFeatureTransformExact(
2.334 + src.first, src.second, src.third,
2.335 + dest.first, dest.second, featurelessPixel);
2.336 +
2.337 +};
2.338 +
2.339 } // namespace enblend
2.340
2.341 #endif /* __NEAREST_H__ */