#include #include #include #include #include #include "TRW-S.h" #define private public #include "typeTruncatedQuadratic2D.h" #undef private #define m_D(pix,l) m_D[(pix)*m_nLabels+(l)] #define m_V(l1,l2) m_V[(l1)*m_nLabels+(l2)] #define MIN(a,b) (((a) < (b)) ? (a) : (b)) #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #define TRUNCATE_MIN(a,b) { if ((a) > (b)) (a) = (b); } #define TRUNCATE_MAX(a,b) { if ((a) < (b)) (a) = (b); } #define TRUNCATE TRUNCATE_MIN ///////////////////////////////////////////////////////////////////////////// // Operations on vectors (arrays of size K) // ///////////////////////////////////////////////////////////////////////////// inline void CopyVector(TRWS::REAL* to, MRF::CostVal* from, int K) { TRWS::REAL* to_finish = to + K; do { *to ++ = *from ++; } while (to < to_finish); } inline void AddVector(TRWS::REAL* to, TRWS::REAL* from, int K) { TRWS::REAL* to_finish = to + K; do { *to ++ += *from ++; } while (to < to_finish); } inline TRWS::REAL SubtractMin(TRWS::REAL *D, int K) { int k; TRWS::REAL delta; delta = D[0]; for (k=1; k=0; k--) { M[k] -= delta; TRUNCATE(M[k], M[k+1] + lambda); TRUNCATE(M[k], lambda*smoothMax); } return delta; } //////////////////////////////////////// // L2 // //////////////////////////////////////// inline TRWS::REAL UpdateMessageL2(TRWS::REAL* M, TRWS::REAL* Di_hat, int K, TRWS::REAL gamma, MRF::CostVal lambda, MRF::CostVal smoothMax, void *buf) { TRWS::REAL* Di = (TRWS::REAL*) buf; int* parabolas = (int*) ((char*)buf + K*sizeof(TRWS::REAL)); int* intersections = parabolas + K; TypeTruncatedQuadratic2D::Edge* tmp = NULL; int k; TRWS::REAL delta; assert(lambda >= 0); Di[0] = gamma*Di_hat[0] - M[0]; delta = Di[0]; for (k=1; kDistanceTransformL2(K, 1, lambda, Di, M, parabolas, intersections); for (k=0; k 0) { for (kj=0; kj= 0); m_type = (smoothExp == 1) ? L1 : L2; int ki, kj; CostVal cost; m_needToFreeV = true; m_V = new CostVal[m_nLabels*m_nLabels]; for (ki=0; ki smoothMax) cost = smoothMax; m_V[ki*m_nLabels + kj] = m_V[kj*m_nLabels + ki] = cost*lambda; } m_smoothMax = smoothMax; m_lambda = lambda; } void TRWS::setCues(CostVal* hCue, CostVal* vCue) { m_horzWeights = hCue; m_vertWeights = vCue; } void TRWS::initializeAlg() { assert(m_type != NONE); int i; // determine type if (m_type == L1 && m_nLabels == 2) { m_type = BINARY; } // allocate messages int messageNum = (m_type == BINARY) ? 4*m_nPixels : 4*m_nPixels*m_nLabels; m_messageArraySizeInBytes = messageNum*sizeof(REAL); m_messages = new REAL[messageNum]; memset(m_messages, 0, messageNum*sizeof(REAL)); if (m_type == BINARY) { assert(m_DBinary == NULL && m_horzWeightsBinary == NULL && m_horzWeightsBinary == NULL); m_DBinary = new CostVal[m_nPixels]; m_horzWeightsBinary = new CostVal[m_nPixels]; m_vertWeightsBinary = new CostVal[m_nPixels]; if ( m_dataType == ARRAY) { for (i=0; i 0) { kj = m_answer[n-1]; for (ki=0; ki 0) { kj = m_answer[n-m_width]; for (ki=0; ki 0) { kj = m_answer[n-1]; for (ki=0; ki 0) { kj = m_answer[n-m_width]; for (ki=0; ki 0) { kj = m_answer[n-1]; CostVal lambda = (m_varWeights) ? m_horzWeights[n-1] : 1; for (ki=0; ki 0) { kj = m_answer[n-m_width]; CostVal lambda = (m_varWeights) ? m_vertWeights[n-m_width] : 1; for (ki=0; ki(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) // compute min delta = Di[0]; m_answer[n] = 0; for (ki=1; ki Di[ki]) { delta = Di[ki]; m_answer[n] = ki; } } } delete [] Di; } else // m_type == BINARY { int x, y, n; REAL* M_ptr; REAL Di; n = 0; M_ptr = m_messages; for (y=0; y 0) Di += (m_answer[n-1] == 0) ? m_horzWeightsBinary[n-1] : -m_horzWeightsBinary[n-1]; if (y > 0) Di += (m_answer[n-m_width] == 0) ? m_vertWeightsBinary[n-m_width] : -m_vertWeightsBinary[n-m_width]; if (x < m_width-1) Di += M_ptr[0]; // message (x+1,y)->(x,y) if (y < m_height-1) Di += M_ptr[1]; // message (x,y+1)->(x,y) // compute min m_answer[n] = (Di >= 0) ? 0 : 1; } } } void TRWS::optimize_GRID_L1(int nIterations) { int x, y, n, K = m_nLabels; CostVal* D_ptr; REAL* M_ptr; REAL* Di; Di = new REAL[K]; for ( ; nIterations > 0; nIterations --) { // forward pass n = 0; D_ptr = m_D; M_ptr = m_messages; for (y=0; y 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) if (x < m_width-1) { CostVal lambda = (m_varWeights) ? m_lambda*m_horzWeights[n] : m_lambda; UpdateMessageL1(M_ptr, Di, K, 0.5, lambda, m_smoothMax); } if (y < m_height-1) { CostVal lambda = (m_varWeights) ? m_lambda*m_vertWeights[n] : m_lambda; UpdateMessageL1(M_ptr+K, Di, K, 0.5, lambda, m_smoothMax); } } // backward pass m_lowerBound = 0; n --; D_ptr -= K; M_ptr -= 2*K; for (y=m_height-1; y>=0; y--) for (x=m_width-1; x>=0; x--, D_ptr-=K, M_ptr-=2*K, n--) { CopyVector(Di, D_ptr, K); if (x > 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) m_lowerBound += SubtractMin(Di, K); if (x > 0) { CostVal lambda = (m_varWeights) ? m_lambda*m_horzWeights[n-1] : m_lambda; m_lowerBound += UpdateMessageL1(M_ptr-2*K, Di, K, 0.5, lambda, m_smoothMax); } if (y > 0) { CostVal lambda = (m_varWeights) ? m_lambda*m_vertWeights[n-m_width] : m_lambda; m_lowerBound += UpdateMessageL1(M_ptr-(2*m_width-1)*K, Di, K, 0.5, lambda, m_smoothMax); } } } delete [] Di; } void TRWS::optimize_GRID_L2(int nIterations) { int x, y, n, K = m_nLabels; CostVal* D_ptr; REAL* M_ptr; REAL* Di; void* buf; Di = new REAL[K]; buf = new char[(2*K+1)*sizeof(int) + K*sizeof(REAL)]; for ( ; nIterations > 0; nIterations --) { // forward pass n = 0; D_ptr = m_D; M_ptr = m_messages; for (y=0; y 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) if (x < m_width-1) { CostVal lambda = (m_varWeights) ? m_lambda*m_horzWeights[n] : m_lambda; UpdateMessageL2(M_ptr, Di, K, 0.5, lambda, m_smoothMax, buf); } if (y < m_height-1) { CostVal lambda = (m_varWeights) ? m_lambda*m_vertWeights[n] : m_lambda; UpdateMessageL2(M_ptr+K, Di, K, 0.5, lambda, m_smoothMax, buf); } } // backward pass m_lowerBound = 0; n --; D_ptr -= K; M_ptr -= 2*K; for (y=m_height-1; y>=0; y--) for (x=m_width-1; x>=0; x--, D_ptr-=K, M_ptr-=2*K, n--) { CopyVector(Di, D_ptr, K); if (x > 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) m_lowerBound += SubtractMin(Di, K); if (x > 0) { CostVal lambda = (m_varWeights) ? m_lambda*m_horzWeights[n-1] : m_lambda; m_lowerBound += UpdateMessageL2(M_ptr-2*K, Di, K, 0.5, lambda, m_smoothMax, buf); } if (y > 0) { CostVal lambda = (m_varWeights) ? m_lambda*m_vertWeights[n-m_width] : m_lambda; m_lowerBound += UpdateMessageL2(M_ptr-(2*m_width-1)*K, Di, K, 0.5, lambda, m_smoothMax, buf); } } } delete [] Di; delete [] (char *)buf; } void TRWS::optimize_GRID_BINARY(int nIterations) { int x, y, n; REAL* M_ptr; REAL Di; for ( ; nIterations > 0; nIterations --) { // forward pass n = 0; M_ptr = m_messages; for (y=0; y 0) Di += M_ptr[-2]; // message (x-1,y)->(x,y) if (y > 0) Di += M_ptr[-2*m_width+1]; // message (x,y-1)->(x,y) if (x < m_width-1) Di += M_ptr[0]; // message (x+1,y)->(x,y) if (y < m_height-1) Di += M_ptr[1]; // message (x,y+1)->(x,y) REAL DiScaled = Di * 0.5; if (x < m_width-1) { Di = DiScaled - M_ptr[0]; CostVal lambda = m_horzWeightsBinary[n]; if (lambda < 0) { Di = -Di; lambda = -lambda; } if (Di > lambda) M_ptr[0] = lambda; else M_ptr[0] = (Di < -lambda) ? -lambda : Di; } if (y < m_height-1) { Di = DiScaled - M_ptr[1]; CostVal lambda = m_vertWeightsBinary[n]; if (lambda < 0) { Di = -Di; lambda = -lambda; } if (Di > lambda) M_ptr[1] = lambda; else M_ptr[1] = (Di < -lambda) ? -lambda : Di; } } // backward pass n --; M_ptr -= 2; for (y=m_height-1; y>=0; y--) for (x=m_width-1; x>=0; x--, M_ptr-=2, n--) { Di = m_DBinary[n]; if (x > 0) Di += M_ptr[-2]; // message (x-1,y)->(x,y) if (y > 0) Di += M_ptr[-2*m_width+1]; // message (x,y-1)->(x,y) if (x < m_width-1) Di += M_ptr[0]; // message (x+1,y)->(x,y) if (y < m_height-1) Di += M_ptr[1]; // message (x,y+1)->(x,y) REAL DiScaled = Di * 0.5; if (x > 0) { Di = DiScaled - M_ptr[-2]; CostVal lambda = m_horzWeightsBinary[n-1]; if (lambda < 0) { Di = -Di; lambda = -lambda; } if (Di > lambda) M_ptr[-2] = lambda; else M_ptr[-2] = (Di < -lambda) ? -lambda : Di; } if (y > 0) { Di = DiScaled - M_ptr[-2*m_width+1]; CostVal lambda = m_vertWeightsBinary[n-m_width]; if (lambda < 0) { Di = -Di; lambda = -lambda; } if (Di > lambda) M_ptr[-2*m_width+1] = lambda; else M_ptr[-2*m_width+1] = (Di < -lambda) ? -lambda : Di; } } } m_lowerBound = 0; } void TRWS::optimize_GRID_FIXED_MATRIX(int nIterations) { int x, y, n, K = m_nLabels; CostVal* D_ptr; REAL* M_ptr; REAL* Di; void* buf; Di = new REAL[K]; buf = new REAL[K]; for ( ; nIterations > 0; nIterations --) { // forward pass n = 0; D_ptr = m_D; M_ptr = m_messages; for (y=0; y 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) if (x < m_width-1) { CostVal lambda = (m_varWeights) ? m_horzWeights[n] : 1; UpdateMessageFIXED_MATRIX(M_ptr, Di, K, 0.5, lambda, m_V, buf); } if (y < m_height-1) { CostVal lambda = (m_varWeights) ? m_vertWeights[n] : 1; UpdateMessageFIXED_MATRIX(M_ptr+K, Di, K, 0.5, lambda, m_V, buf); } } // backward pass m_lowerBound = 0; n --; D_ptr -= K; M_ptr -= 2*K; for (y=m_height-1; y>=0; y--) for (x=m_width-1; x>=0; x--, D_ptr-=K, M_ptr-=2*K, n--) { CopyVector(Di, D_ptr, K); if (x > 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) m_lowerBound += SubtractMin(Di, K); if (x > 0) { CostVal lambda = (m_varWeights) ? m_horzWeights[n-1] : 1; m_lowerBound += UpdateMessageFIXED_MATRIX(M_ptr-2*K, Di, K, 0.5, lambda, m_V, buf); } if (y > 0) { CostVal lambda = (m_varWeights) ? m_vertWeights[n-m_width] : 1; m_lowerBound += UpdateMessageFIXED_MATRIX(M_ptr-(2*m_width-1)*K, Di, K, 0.5, lambda, m_V, buf); } } } delete [] Di; delete [] (REAL *)buf; } void TRWS::optimize_GRID_GENERAL(int nIterations) { int x, y, n, K = m_nLabels; CostVal* D_ptr; REAL* M_ptr; REAL* Di; void* buf; Di = new REAL[K]; buf = new REAL[K]; for ( ; nIterations > 0; nIterations --) { // forward pass n = 0; D_ptr = m_D; M_ptr = m_messages; CostVal* V_ptr = m_V; for (y=0; y 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) if (x < m_width-1) { if (m_V) UpdateMessageGENERAL(M_ptr, Di, K, 0.5, /* forward dir*/ 0, V_ptr, buf); else UpdateMessageGENERAL(M_ptr, Di, K, 0.5, m_smoothFn, n, n+1, buf); } if (y < m_height-1) { if (m_V) UpdateMessageGENERAL(M_ptr+K, Di, K, 0.5, /* forward dir*/ 0, V_ptr+K*K, buf); else UpdateMessageGENERAL(M_ptr+K, Di, K, 0.5, m_smoothFn, n, n+m_width, buf); } } // backward pass m_lowerBound = 0; n --; D_ptr -= K; M_ptr -= 2*K; V_ptr -= 2*K*K; for (y=m_height-1; y>=0; y--) for (x=m_width-1; x>=0; x--, D_ptr-=K, M_ptr-=2*K, V_ptr-=2*K*K, n--) { CopyVector(Di, D_ptr, K); if (x > 0) AddVector(Di, M_ptr-2*K, K); // message (x-1,y)->(x,y) if (y > 0) AddVector(Di, M_ptr-(2*m_width-1)*K, K); // message (x,y-1)->(x,y) if (x < m_width-1) AddVector(Di, M_ptr, K); // message (x+1,y)->(x,y) if (y < m_height-1) AddVector(Di, M_ptr+K, K); // message (x,y+1)->(x,y) // normalize Di, update lower bound m_lowerBound += SubtractMin(Di, K); if (x > 0) { if (m_V) m_lowerBound += UpdateMessageGENERAL(M_ptr-2*K, Di, K, 0.5, /* backward dir */ 1, V_ptr-2*K*K, buf); else m_lowerBound += UpdateMessageGENERAL(M_ptr-2*K, Di, K, 0.5, m_smoothFn, n, n-1, buf); } if (y > 0) { if (m_V) m_lowerBound += UpdateMessageGENERAL(M_ptr-(2*m_width-1)*K, Di, K, 0.5, /* backward dir */ 1, V_ptr-(2*m_width-1)*K*K, buf); else m_lowerBound += UpdateMessageGENERAL(M_ptr-(2*m_width-1)*K, Di, K, 0.5, m_smoothFn, n, n-m_width, buf); } } } delete [] Di; delete [] (REAL *)buf; }