Total Variation Deconvolution using Split Bregman
dsolve_inc.c
Go to the documentation of this file.
1 
15 #include "tvregopt.h"
16 
17 
47 static void DSolve(tvregsolver *S)
48 {
49  numvec2 *d = S->d;
50  numvec2 *dtilde = S->dtilde;
51  const num *u = S->u;
52  const int Width = S->Width;
53  const int Height = S->Height;
54  const int NumChannels = S->NumChannels;
55  const num Thresh = 1/S->Opt.Gamma1;
56  const num ThreshSquared = Thresh * Thresh;
57  const long ChannelStride = ((long)Width) * ((long)Height);
58  const long NumEl = NumChannels * ChannelStride;
59  numvec2 dnew;
60  num Magnitude;
61  long i;
62  int x, y;
63 
64  for(y = 0; y < Height - 1; y++, d++, dtilde++, u++)
65  {
66  /* Perform vectorial shrinkage for interior points */
67  for(x = 0; x < Width - 1; x++, d++, dtilde++, u++)
68  {
69  for(i = 0, Magnitude = 0; i < NumEl; i += ChannelStride)
70  {
71  d[i].x += (u[i + 1] - u[i]) - dtilde[i].x;
72  d[i].y += (u[i + Width] - u[i]) - dtilde[i].y;
73  Magnitude += d[i].x*d[i].x + d[i].y*d[i].y;
74  }
75 
76  if(Magnitude > ThreshSquared)
77  {
78  Magnitude = 1 - Thresh/(num)sqrt(Magnitude);
79 
80  for(i = 0; i < NumEl; i += ChannelStride)
81  {
82  dnew.x = Magnitude*d[i].x;
83  dnew.y = Magnitude*d[i].y;
84  dtilde[i].x = 2*dnew.x - d[i].x;
85  dtilde[i].y = 2*dnew.y - d[i].y;
86  d[i] = dnew;
87  }
88  }
89  else
90  for(i = 0; i < NumEl; i += ChannelStride)
91  {
92  dtilde[i].x = -d[i].x;
93  dtilde[i].y = -d[i].y;
94  d[i].x = 0;
95  d[i].y = 0;
96  }
97  }
98 
99  /* Right edge */
100  for(i = 0, Magnitude = 0; i < NumEl; i += ChannelStride)
101  {
102  d[i].y += (u[i + Width] - u[i]) - dtilde[i].y;
103  Magnitude += d[i].y*d[i].y;
104  d[i].x = dtilde[i].x = 0;
105  }
106 
107  if(Magnitude > ThreshSquared)
108  {
109  Magnitude = 1 - Thresh/(num)sqrt(Magnitude);
110 
111  for(i = 0; i < NumEl; i += ChannelStride)
112  {
113  dnew.y = Magnitude*d[i].y;
114  dtilde[i].y = 2*dnew.y - d[i].y;
115  d[i].y = dnew.y;
116  }
117  }
118  else
119  for(i = 0; i < NumEl; i += ChannelStride)
120  {
121  dtilde[i].y = -d[i].y;
122  d[i].y = 0;
123  }
124  }
125 
126  /* Bottom edge */
127  for(x = 0; x < Width - 1; x++, d++, dtilde++, u++)
128  {
129  for(i = 0, Magnitude = 0; i < NumEl; i += ChannelStride)
130  {
131  d[i].x += (u[i + 1] - u[i]) - dtilde[i].x;
132  Magnitude += d[i].x*d[i].x;
133  d[i].y = dtilde[i].y = 0;
134  }
135 
136  if(Magnitude > ThreshSquared)
137  {
138  Magnitude = 1 - Thresh/(num)sqrt(Magnitude);
139 
140  for(i = 0; i < NumEl; i += ChannelStride)
141  {
142  dnew.x = Magnitude*d[i].x;
143  dtilde[i].x = 2*dnew.x - d[i].x;
144  d[i].x = dnew.x;
145  }
146  }
147  else
148  for(i = 0; i < NumEl; i += ChannelStride)
149  {
150  dtilde[i].x = -d[i].x;
151  d[i].x = 0;
152  }
153  }
154 
155  /* Bottom-right corner */
156  for(i = 0; i < NumEl; i += ChannelStride)
157  d[i].x = d[i].y = dtilde[i].x = dtilde[i].y = 0;
158 }