tensor3d.cpp

把黑白照片转换成彩色照片的数值优化算法。

  for(y = 1 ; y < r.len_y - 1 ; y++) {
    for( x = 1 ; x < r.len_x - 1 ; x++) {
      i = x + y * r.len_x ;
      j = 2 * x + 2 * y * len_x ;

      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-1] + data[j+1] + data[j-len_x] + data[j+len_x] + data[j+plane]) ;

      i += (r.len_z - 1) * r.plane ;
      j += (len_z - 1) * plane ;
      
      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-1] + data[j+1] + data[j-len_x] + data[j+len_x] + data[j-plane]) ;
    }
  }

  /* Y-Plane (ZX) */
  
  for( z = 1 ; z < r.len_z - 1 ; z++) {
    for( x = 1 ; x < r.len_x - 1 ; x++) {
      i = x + z * r.plane ;
      j = 2 * x + 2 * z * plane ;

      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-1] + data[j+1] + data[j-plane] + data[j+plane] + data[j+len_x]) ;

      i += (r.len_y-1) * r.len_x ;
      j += (len_y-1) * len_x ;
      
      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-1] + data[j+1] + data[j-plane] + data[j+plane] + data[j-len_x]) ;
    }
  }

  /* X-Plane (ZY) */
  
  for( z = 1 ; z < r.len_z - 1 ; z++) {
    for( y = 1 ; y < r.len_y - 1 ; y++) {
      i = y * r.len_x + z * r.plane ;
      j = 2 * y * len_x + 2 * z * plane ;

      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-len_x] + data[j+len_x] + data[j-plane] + data[j+plane] + data[j+1]) ;

      i += r.len_x - 1 ;
      j += len_x - 1 ;
      
      r.data[i] = 0.2857 * data[j]
	+ 0.14285 * (data[j-len_x] + data[j+len_x] + data[j-plane] + data[j+plane] + data[j-1]) ;
    }
  }

  /* X-Axis */
  
  for( x = 1 ; x < r.len_x - 1 ; x++) {
    i = x ;
    j = 2 * x ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+1] + data[j-1] + data[j+len_x] + data[j+plane]) ;

    i = x + (r.len_y-1) * r.len_x ;
    j = 2 * x + (len_y-1) * len_x ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+1] + data[j-1] + data[j-len_x] + data[j+plane]) ;
    
    i = x + (r.len_z-1) * r.plane ;
    j = 2 * x + (len_z-1) * plane ;
    
    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+1] + data[j-1] + data[j+len_x] + data[j-plane]) ;

    i = x + (r.len_z-1) * r.plane + (r.len_y-1) * r.len_x ;
    j = 2 * x + (len_z-1) * plane + (len_y-1) * len_x ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+1] + data[j-1] + data[j-len_x] + data[j-plane]) ;

  }

  /* Y-Axis */

  for( y = 1 ; y < r.len_y - 1 ; y++) {
    
    i = y * r.len_x ;
    j = 2 * y * len_x ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+len_x] + data[j-len_x] + data[j+1] + data[j+plane]) ;

    i = y * r.len_x + r.len_x - 1 ;
    j = 2 * y * len_x + len_x - 1 ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+len_x] + data[j-len_x] + data[j-1] + data[j+plane]) ;

    i = y * r.len_x + (r.len_y-1) * r.plane ;
    j = 2 * y * len_x + (len_y-1) * plane ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+len_x] + data[j-len_x] + data[j+1] + data[j-plane]) ;
    
    i = y * r.len_x + (r.len_y-1) * r.plane + r.len_x - 1 ;
    j = 2 * y * len_x + (len_y-1) * plane + len_x - 1 ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j-len_x] + data[j-len_x] + data[j+1] + data[j-plane]) ;
  }
  
  /* Z-Axis */
     
  for( z = 1 ; z < r.len_z - 1 ; z++) {
    i = z * r.plane ;
    j = 2 * z * plane ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+plane] + data[j-plane] + data[j+len_x] + data[j+1]) ;
  
    i = z * r.plane + r.len_x - 1 ;
    j = 2 * z * plane +len_x - 1 ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+plane] + data[j-plane] + data[j+len_x] + data[j-1]) ;

    i = z * r.plane + (r.len_y-1) * r.len_x ;
    j = 2 * z * plane + (len_y-1) * len_x ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+plane] + data[j-plane] + data[j-len_x] + data[j+1]) ;

    i = z * r.plane + (r.len_y-1) * r.len_x + r.len_x - 1 ;
    j = 2 * z * plane + (len_y-1) * len_x + len_x - 1 ;

    r.data[i] = 0.3333 * data[j]
      + 0.16666 * (data[j+plane] + data[j-plane] + data[j-len_x] + data[j-1]) ;
  }

  i = 0 ;
  j = 0 ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j+1] + data[j+len_x] + data[j+plane]) ;

  i = r.len_x - 1 ;
  j = len_x - 1 ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j-1] + data[j+len_x] + data[j+plane]) ;

  i = (r.len_y - 1) * r.len_x ;
  j = (len_y - 1) * len_x  ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j+1] + data[j-len_x] + data[j+plane]) ;
  
  i = (r.len_y - 1) * r.len_x + r.len_x - 1 ; 
  j = (len_y - 1) * len_x + len_x - 1 ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j-1] + data[j-len_x] + data[j+plane]) ;

  i = (r.len_z - 1) * r.plane ;
  j = (len_z - 1 ) * plane ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j+1] + data[j+len_x] + data[j-plane]) ;

  i = r.len_x - 1 + (r.len_z - 1) * r.plane ;
  j = len_x - 1 + (len_z - 1 ) * plane ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j-1] + data[j+len_x] + data[j-plane]) ;

  i = (r.len_y - 1) * r.len_x + (r.len_z - 1) * r.plane ;
  j = (len_y - 1) * len_x + (len_z - 1 ) * plane  ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j+1] + data[j-len_x] + data[j-plane]) ;
  
  i = (r.len_y - 1) * r.len_x + r.len_x - 1 + (r.len_z - 1) * r.plane ; 
  j = (len_y - 1) * len_x + len_x - 1 + (len_z - 1 ) * plane ;
  
  r.data[i] = 0.4 * data[j]
    + 0.2 * (data[j-1] + data[j-len_x] + data[j-plane]) ;
}

void Tensor3d::restrict_none(Tensor3d& r) {
  r.set((len_x - 1) / 2 + 1, (len_y - 1) / 2 + 1, (len_z - 1) / 2 + 1) ;

  int i, j,x,y,z ;
   
  for( z = 0 ; z < r.len_z ; z++) {
    for(y = 0 ; y < r.len_y ; y++) {
      for( x = 0 ; x < r.len_x ; x++) {
	i = x + y * r.len_x + z * r.plane ;
	j = 2 * x + 2 * y * len_x + 2 * z * plane ;

	r.data[i] = data[j] ;
      }
    }
  }
}

void Tensor3d::prolong_full(Tensor3d& p) {
  p.set((len_x - 1) * 2 + 1, (len_y - 1) * 2 + 1, (len_z - 1) * 2 + 1) ;

  int i, j ,x,y,z;
  
  p.data[0] = data[0] ;
  
  for(x = 1 ; x < len_x ; x++) {
    i = x ;
    j = 2 * x ;

    p.data[j] = data[i] ;
    p.data[j-1] = 0.5 * (data[i] + data[i-1]) ;
  }
  
  for( y = 1 ; y < len_y ; y++) {
    i = y * len_x ;
    j = 2 * y * p.len_x  ;

    p.data[j] = data[i] ;
    p.data[j-p.len_x] = 0.5 * (data[i] + data[i-len_x]) ;
  }

  
  for( z = 1 ; z < len_z ; z++) {
    i = z * plane ;
    j = 2 * z * p.plane  ;

    p.data[j] = data[i] ;
    p.data[j-p.plane] = 0.5 * (data[i] + data[i-plane]) ;
  }

  for( y = 1 ; y < len_y ; y++) {
    for( x = 1 ; x < len_x ; x++) {
      i = x + y * len_x ;
      j = 2 * x + 2 * y * p.len_x  ;
      
      p.data[j] = data[i] ;
      
      p.data[j-1] = 0.5 * (data[i] + data[i-1]) ;
      p.data[j-p.len_x] = 0.5 * (data[i] + data[i-len_x]) ;

      p.data[j-1-p.len_x] = 0.25 * (data[i] + data[i-len_x] + data[i-1] + data[i-1-len_x]) ;
    }
  }
  
  for(z = 1 ; z < len_z ; z++) {
    for( x = 1 ; x < len_x ; x++) {
      i = x + z * plane ;
      j = 2 * x + 2 * z * p.plane  ;
      
      p.data[j] = data[i] ;
      
      p.data[j-1] = 0.5 * (data[i] + data[i-1]) ;
      p.data[j-p.plane] = 0.5 * (data[i] + data[i-plane]) ;

      p.data[j-1-plane] = 0.25 * (data[i] + data[i-plane] + data[i-1] + data[i-1-plane]) ;
    }
  }
  
  for( z = 1 ; z < len_z ; z++) {
    for( y = 1 ; y < len_y ; y++) {
      i = y * len_x + z * plane ;
      j = 2 * y * p.len_x + 2 * z * p.plane  ;
      
      p.data[j] = data[i] ;
      
      p.data[j-p.len_x] = 0.5 * (data[i] + data[i-len_x]) ;
      p.data[j-p.plane] = 0.5 * (data[i] + data[i-plane]) ;

      p.data[j-p.len_x-plane] = 0.25 * (data[i] + data[i-plane] + data[i-len_x] + data[i-len_x-plane]) ;
    }
  }
    
  for( z = 1 ; z < len_z ; z++) {
    for( y = 1 ; y < len_y ; y++) {
      for( x = 1 ; x < len_x ; x++) {
	i = x + y * len_x + z * plane ;
	j = 2 * x + 2 * y * p.len_x + 2 * z * p.plane ;

	p.data[j] = data[i] ;

	p.data[j-1] = 0.5 * (data[i-1] + data[i]) ;
	p.data[j-p.len_x] = 0.5 * (data[i-len_x] + data[i]) ;
	p.data[j-p.plane] = 0.5 * (data[i-plane] + data[i]) ;

	p.data[j-1-p.len_x] = 0.25 * (data[i-1] + data[i] + data[i-len_x] + data[i-len_x-1]) ;
	p.data[j-1-p.plane] = 0.25 * (data[i-1] + data[i] + data[i-plane] + data[i-plane-1]) ;
	p.data[j-p.len_x-p.plane] = 0.25 * (data[i-len_x] + data[i] + data[i-plane] + data[i-plane-len_x]) ;

	p.data[j-1-p.len_x-p.plane] = 0.125 * (data[i-len_x] + data[i] + data[i-plane] + data[i-plane-len_x] +
					       data[i-len_x-1] + data[i-1-plane] + data[i-1] + data[i-plane-len_x-1]) ;
      }
    }
  }
}

void Tensor3d::prolong_none(Tensor3d& p) {
  p.set((len_x - 1) * 2 + 1, (len_y - 1) * 2 + 1, (len_z - 1) * 2 + 1) ;

  int i, j,x,y,z ;
  
  p.data[0] = data[0] ;
  
  for(x = 1 ; x < len_x ; x++) {
    i = x ;
    j = 2 * x ;

    p.data[j-1] = p.data[j] = data[i] ;
  }
  
  for(y = 1 ; y < len_y ; y++) {
    i = y * len_x ;
    j = 2 * y * p.len_x  ;

    p.data[j-p.len_x] = p.data[j] = data[i] ;
  }

  
  for( z = 1 ; z < len_z ; z++) {
    i = z * plane ;
    j = 2 * z * p.plane  ;

    p.data[j-p.plane] = p.data[j] = data[i] ;
  }

  for( y = 1 ; y < len_y ; y++) {
    for( x = 1 ; x < len_x ; x++) {
      i = x + y * len_x ;
      j = 2 * x + 2 * y * p.len_x  ;
      
      p.data[j-1-p.len_x] = p.data[j-p.len_x] = p.data[j-1] = p.data[j] = data[i] ;
    }
  }
  
  for( z = 1 ; z < len_z ; z++) {
    for( x = 1 ; x < len_x ; x++) {
      i = x + z * plane ;
      j = 2 * x + 2 * z * p.plane  ;
      
      p.data[j-1-plane] = p.data[j-p.plane] = p.data[j-1] = p.data[j] = data[i] ;
    }
  }
  
  for( z = 1 ; z < len_z ; z++) {
    for( y = 1 ; y < len_y ; y++) {
      i = y * len_x + z * plane ;
      j = 2 * y * p.len_x + 2 * z * p.plane  ;
      
      p.data[j-p.len_x-plane] = p.data[j-p.plane] = p.data[j-p.len_x] = p.data[j] = data[i] ;
    }
  }
    
  for( z = 1 ; z < len_z ; z++) {
    for( y = 1 ; y < len_y ; y++) {
      for( x = 1 ; x < len_x ; x++) {
	i = x + y * len_x + z * plane ;
	j = 2 * x + 2 * y * p.len_x + 2 * z * p.plane ;
	
	p.data[j-1-p.len_x-p.plane] = p.data[j-p.len_x-p.plane] = p.data[j-1-p.plane] = p.data[j-1-p.len_x] = p.data[j-p.plane] = p.data[j-p.len_x] = p.data[j-1] = p.data[j] = data[i] ;
      }
    }
  }
}

  
double Tensor3d::norm_max() {
  double norm = 0 ;

  for(int i = 0 ; i < length ; i++)
    if(fabs(data[i]) > norm)
      norm = fabs(data[i]) ;

  return norm ;
}

void Tensor3d::one_minus() {
  for(int i = 0 ; i < length ; i++)
    data[i] = 1.0 - data[i] ;
}
 
Tensor3d& Tensor3d::operator*=(real val) {
  for(int i = 0 ; i < length ; i++)
     data[i] *= val ;
  return *this ;
}

Tensor3d& Tensor3d::operator+=(Tensor3d& other) {
  if(length != other.length) {
    cerr << "Vectors length mismatch! (op+)" << endl ;
    exit(1) ;
  }
  
  for(int i = 0 ; i < length ; i++)
    data[i] += other.data[i] ;
  return *this ;  
}

Tensor3d& Tensor3d::operator-=(Tensor3d& other) {
  if(length != other.length) {
    cerr << "Vectors length mismatch! (op-)" << endl ;
    exit(1) ;
  }
  
  for(int i = 0 ; i < length ; i++)
    data[i] -= other.data[i] ;
  return *this ;
}

double Tensor3d::norm_abs() {
  double norm = 0 ;
  
  for (int i = 0 ; i < length ; i++)
    norm += fabs(data[i]) ;

  return norm / length ;
}

double Tensor3d::norm_euc() {
  double norm = 0 ;
  
  for (int i = 0 ; i < length ; i++)
    norm += data[i] * data[i] ;

  //return norm / length ;
  return norm ;
}