Hi Francois,<br><br>> that is to say either I1 or I2 have size 1 on one of<br>
> the coordinates. <br><br>I'm not sure that should be the case when you are<br>computing the gradient for a 2D image.<br><br>> Assuming I do a 2nd derivative, I should request a region<br>
> padded by a radius of 2 right ?<br><br>Yes, assuming that you are using central differences<br>to compute the second derivative, you will need at<br>least 5 values (which is a radius = 2).<br><br><br> Regards,<br>
<br><br> Luis<br><br><br>-------------------------------------------------------------------------<br><div class="gmail_quote">On Thu, Sep 23, 2010 at 12:13 PM, <span dir="ltr"><<a href="mailto:devieill@irit.fr">devieill@irit.fr</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">Hi all,<br>
<br>
I am currently rewriting some of my code, and as a matter of fact it find<br>
it iteresting to write a filter to easily compute the gradient of the<br>
total variation. I have already done this computation with a standard itk<br>
pipeline, at a high cost (swapping) for large image.<br>
So starting from itkGradientMagnitudeImageFilter I have decided to give it<br>
a try. Right now I am bit stuck at computing the properly the second<br>
derivative pass. To explain better from an image I with 2 dimensions, one<br>
has to get the partial derivatives dI/dx, dI/dy and then to form 2 images<br>
:<br>
I1 = (dI/dx)/( (dI/dx)^2 + (dI/dy)^2 + beta^2)<br>
I2 = (dI/dy)/( (dI/dx)^2 + (dI/dy)^2 + beta^2)<br>
And finally to compute the result R:<br>
R = dI1/dx + dI2/dy<br>
<br>
this part is a bit tricky to me since I am not sure to have the proper<br>
images to compute it, that is to say either I1 or I2 have size 1 on one of<br>
the coordinates. Assuming I do a 2nd derivative, I should request a region<br>
padded by a radius of 2 right ?<br>
<br>
However I am not too sure about how to doing it properly and more<br>
importantly if that is the critical point. Here is the important part of<br>
the code doing this.<br>
<br>
Any help would be much appreciated !!<br>
<br>
Regards,<br>
de Vieilleville François<br>
<br>
<br>
template< class TInputImage, class TOutputImage><br>
void<br>
GradientTotalVariationL2ApproximationImageFilter< TInputImage, TOutputImage ><br>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,<br>
int threadId)<br>
{<br>
//////////////////////////////////////////////<br>
// create copies for storing the returned values<br>
typename TOutputImage::Pointer ddx[ImageDimension];<br>
<br>
//allocate copies<br>
for (unsigned int i = 0; i < ImageDimension; ++i) {<br>
ddx[i] = TOutputImage::New();<br>
ddx[i]->SetRegions(outputRegionForThread );<br>
ddx[i]->Allocate();<br>
}<br>
<br>
OutputPixelType grad_tv;<br>
ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;<br>
<br>
ConstNeighborhoodIterator<InputImageType> nit;<br>
ConstNeighborhoodIterator<TInputImage> bit;<br>
ImageRegionIterator<OutputImageType> it;<br>
<br>
NeighborhoodInnerProduct<InputImageType, RealType> SIP;<br>
<br>
// Get the input and output<br>
OutputImageType * outputImage = this->GetOutput();<br>
const InputImageType * inputImage = this->GetInput();<br>
<br>
// Set up operators<br>
DerivativeOperator<RealType, ImageDimension> op[ImageDimension];<br>
<br>
for (int i = 0; i< ImageDimension; i++)<br>
{<br>
op[i].SetDirection(0);<br>
op[i].SetOrder(1);<br>
op[i].CreateDirectional();<br>
<br>
// Reverse order of coefficients for the convolution with the image to<br>
// follow.<br>
//op[i].FlipAxes();<br>
<br>
// Take into account the pixel spacing if necessary<br>
if (m_UseImageSpacing == true)<br>
{<br>
if ( this->GetInput()->GetSpacing()[i] == 0.0 )<br>
{<br>
itkExceptionMacro(<< "Image spacing cannot be zero.");<br>
}<br>
else<br>
{<br>
op[i].ScaleCoefficients( 1.0 / this->GetInput()->GetSpacing()[i] );<br>
}<br>
}<br>
}<br>
<br>
// Calculate iterator radius<br>
Size<ImageDimension> radius;<br>
for (unsigned int i = 0; i < ImageDimension; ++i)<br>
{<br>
radius[i] = op[0].GetRadius()[0];<br>
//std::cout << "radius["<<i<<"]=" << radius[i] << std::endl;<br>
}<br>
<br>
// Find the data-set boundary "faces"<br>
typename<br>
NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType<br>
faceList;<br>
NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage> bC;<br>
faceList = bC(inputImage, outputRegionForThread, radius);<br>
<br>
typename<br>
NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator<br>
fit;<br>
fit = faceList.begin();<br>
<br>
// support progress methods/callbacks<br>
ProgressReporter progress(this, threadId,<br>
outputRegionForThread.GetNumberOfPixels());<br>
<br>
// Initialize the x_slice array<br>
nit = ConstNeighborhoodIterator<InputImageType>(radius, inputImage, *fit);<br>
<br>
//////////////////////////////////////////////<br>
<br>
std::slice x_slice[ImageDimension];<br>
const unsigned long center = nit.Size() / 2;<br>
for (unsigned int i = 0; i < ImageDimension; ++i)<br>
{<br>
<br>
x_slice[i] = std::slice( center - nit.GetStride(i) * radius[0],<br>
op[i].GetSize()[0], nit.GetStride(i));<br>
}<br>
// Process non-boundary face and then each of the boundary faces.<br>
// These are N-d regions which border the edge of the buffer.<br>
<br>
// set the proper iterators on the copied regions<br>
ImageRegionIterator<OutputImageType> it_copy[ImageDimension];<br>
// typename TInputImage::RegionType region_on_copy_current_face;<br>
<br>
for (fit = faceList.begin(); fit != faceList.end(); ++fit)<br>
{<br>
bit = ConstNeighborhoodIterator<InputImageType>(radius,<br>
inputImage, *fit);<br>
it = ImageRegionIterator<OutputImageType>(outputImage, *fit);<br>
// set the copied regions iterators<br>
for (unsigned int j =0; j < ImageDimension; ++j) {<br>
it_copy[j] = ImageRegionIterator<InputImageType>(ddx[j], *fit);<br>
}<br>
<br>
bit.OverrideBoundaryCondition(&nbc);<br>
bit.GoToBegin();<br>
<br>
while ( ! bit.IsAtEnd() )<br>
{<br>
RealType a = NumericTraits<RealType>::Zero;<br>
RealType deriv[ImageDimension];<br>
for (int i = 0; i < ImageDimension; ++i)<br>
{<br>
const RealType g = SIP(x_slice[i], bit, op[i]);<br>
deriv[i] = g;<br>
a += g * g;<br>
// vcl_pow(val,2.);<br>
}<br>
it.Value() = vcl_sqrt(a + m_BetaSquare);<br>
for (unsigned int j =0; j < ImageDimension; ++j) {<br>
it_copy[j].Value() = deriv[j]/it.Get();<br>
++(it_copy[j]);<br>
}<br>
++bit;<br>
++it;<br>
}<br>
}<br>
<br>
faceList = bC(ddx[0], outputRegionForThread, radius); // assuming they<br>
are all the same for the ddx[i]....<br>
fit = faceList.begin();<br>
ConstNeighborhoodIterator<InputImageType> nit_copy[ImageDimension];<br>
ConstNeighborhoodIterator<InputImageType> bit_copy[ImageDimension];<br>
<br>
for ( unsigned int i =0; i < ImageDimension; ++i) {<br>
nit_copy[i] = ConstNeighborhoodIterator<InputImageType>(radius,<br>
ddx[i], *fit);<br>
}<br>
// Initialize the x_slice_copy array<br>
std::slice x_slice_copy[ImageDimension];<br>
for (unsigned int i = 0; i < ImageDimension; ++i) {<br>
const unsigned long center = nit_copy[i].Size() / 2;<br>
x_slice_copy[i] = std::slice( center - nit_copy[i].GetStride(i) *<br>
radius[0],<br>
op[i].GetSize()[0], nit_copy[i].GetStride(i));<br>
}<br>
<br>
for (fit = faceList.begin(); fit != faceList.end(); ++fit)<br>
{<br>
for ( unsigned int i =0; i < ImageDimension; ++i) {<br>
bit_copy[i] = ConstNeighborhoodIterator<InputImageType>(radius,<br>
ddx[i], *fit);<br>
bit_copy[i].OverrideBoundaryCondition(&nbc);<br>
bit_copy[i].GoToBegin();<br>
}<br>
it = ImageRegionIterator<OutputImageType>(outputImage, *fit);<br>
//bool nit_copy_are_at_end = false;<br>
while ( !bit_copy[0].IsAtEnd() )<br>
{<br>
grad_tv =0.;<br>
for (unsigned int i = 0; i < ImageDimension; ++i)<br>
{<br>
// compute derivative of previous computed region<br>
double val = SIP(x_slice_copy[i], bit_copy[i], op[i]);<br>
grad_tv += val;<br>
}<br>
it.Value() = grad_tv;<br>
for (unsigned int j =0; j < ImageDimension; ++j) {<br>
++(bit_copy[j]);<br>
}<br>
++it;<br>
progress.CompletedPixel();<br>
}<br>
}<br>
}<br>
<br>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
</blockquote></div><br>