Mercurial > hg > orthanc
comparison Core/Images/ImageProcessing.cpp @ 3503:46cf170ba121
ImageProcessing::SeparableConvolution()
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 27 Aug 2019 12:10:13 +0200 |
parents | c160eafc42a9 |
children | 18566f9e1831 |
comparison
equal
deleted
inserted
replaced
3502:c160eafc42a9 | 3503:46cf170ba121 |
---|---|
1763 | 1763 |
1764 default: | 1764 default: |
1765 throw OrthancException(ErrorCode_NotImplemented); | 1765 throw OrthancException(ErrorCode_NotImplemented); |
1766 } | 1766 } |
1767 } | 1767 } |
1768 | |
1769 | |
1770 // This is a slow implementation of horizontal convolution on one | |
1771 // individual channel, that checks for out-of-image values | |
1772 template <typename RawPixel, unsigned int ChannelsCount> | |
1773 static float GetHorizontalConvolutionFloatSecure(const Orthanc::ImageAccessor& source, | |
1774 const std::vector<float>& horizontal, | |
1775 size_t horizontalAnchor, | |
1776 unsigned int x, | |
1777 unsigned int y, | |
1778 float leftBorder, | |
1779 float rightBorder, | |
1780 unsigned int channel) | |
1781 { | |
1782 const RawPixel* row = reinterpret_cast<const RawPixel*>(source.GetConstRow(y)) + channel; | |
1783 | |
1784 float p = 0; | |
1785 | |
1786 for (unsigned int k = 0; k < horizontal.size(); k++) | |
1787 { | |
1788 float value; | |
1789 | |
1790 if (x + k < horizontalAnchor) // Negation of "x - horizontalAnchor + k >= 0" | |
1791 { | |
1792 value = leftBorder; | |
1793 } | |
1794 else if (x + k >= source.GetWidth() + horizontalAnchor) // Negation of "x - horizontalAnchor + k < width" | |
1795 { | |
1796 value = rightBorder; | |
1797 } | |
1798 else | |
1799 { | |
1800 // The value lies within the image | |
1801 value = row[(x - horizontalAnchor + k) * ChannelsCount]; | |
1802 } | |
1803 | |
1804 p += value * horizontal[k]; | |
1805 } | |
1806 | |
1807 return p; | |
1808 } | |
1809 | |
1810 | |
1811 | |
1812 // This is an implementation of separable convolution that uses | |
1813 // floating-point arithmetics, and an intermediate Float32 | |
1814 // image. The out-of-image values are taken as the border | |
1815 // value. Further optimization is possible. | |
1816 template <typename RawPixel, unsigned int ChannelsCount> | |
1817 static void SeparableConvolutionFloat(ImageAccessor& image /* inplace */, | |
1818 const std::vector<float>& horizontal, | |
1819 size_t horizontalAnchor, | |
1820 const std::vector<float>& vertical, | |
1821 size_t verticalAnchor, | |
1822 float normalization) | |
1823 { | |
1824 const unsigned int width = image.GetWidth(); | |
1825 const unsigned int height = image.GetHeight(); | |
1826 | |
1827 | |
1828 /** | |
1829 * Horizontal convolution | |
1830 **/ | |
1831 | |
1832 Image tmp(PixelFormat_Float32, ChannelsCount * width, height, false); | |
1833 | |
1834 for (unsigned int y = 0; y < height; y++) | |
1835 { | |
1836 const RawPixel* row = reinterpret_cast<const RawPixel*>(image.GetConstRow(y)); | |
1837 | |
1838 float leftBorder[ChannelsCount], rightBorder[ChannelsCount]; | |
1839 | |
1840 for (unsigned int c = 0; c < ChannelsCount; c++) | |
1841 { | |
1842 leftBorder[c] = row[c]; | |
1843 rightBorder[c] = row[ChannelsCount * (width - 1) + c]; | |
1844 } | |
1845 | |
1846 float* p = static_cast<float*>(tmp.GetRow(y)); | |
1847 | |
1848 if (width < horizontal.size()) | |
1849 { | |
1850 // It is not possible to have the full kernel within the image, use the direct implementation | |
1851 for (unsigned int x = 0; x < width; x++) | |
1852 { | |
1853 for (unsigned int c = 0; c < ChannelsCount; c++, p++) | |
1854 { | |
1855 *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount> | |
1856 (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c); | |
1857 } | |
1858 } | |
1859 } | |
1860 else | |
1861 { | |
1862 // Deal with the left border | |
1863 for (unsigned int x = 0; x < horizontalAnchor; x++) | |
1864 { | |
1865 for (unsigned int c = 0; c < ChannelsCount; c++, p++) | |
1866 { | |
1867 *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount> | |
1868 (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c); | |
1869 } | |
1870 } | |
1871 | |
1872 // Deal with the central portion of the image (all pixel values | |
1873 // scanned by the kernel lie inside the image) | |
1874 | |
1875 for (unsigned int x = 0; x < width - horizontal.size() + 1; x++) | |
1876 { | |
1877 for (unsigned int c = 0; c < ChannelsCount; c++, p++) | |
1878 { | |
1879 *p = 0; | |
1880 for (unsigned int k = 0; k < horizontal.size(); k++) | |
1881 { | |
1882 *p += static_cast<float>(row[(x + k) * ChannelsCount + c]) * horizontal[k]; | |
1883 } | |
1884 } | |
1885 } | |
1886 | |
1887 // Deal with the right border | |
1888 for (unsigned int x = horizontalAnchor + width - horizontal.size() + 1; x < width; x++) | |
1889 { | |
1890 for (unsigned int c = 0; c < ChannelsCount; c++, p++) | |
1891 { | |
1892 *p = GetHorizontalConvolutionFloatSecure<RawPixel, ChannelsCount> | |
1893 (image, horizontal, horizontalAnchor, x, y, leftBorder[c], rightBorder[c], c); | |
1894 } | |
1895 } | |
1896 } | |
1897 } | |
1898 | |
1899 | |
1900 /** | |
1901 * Vertical convolution | |
1902 **/ | |
1903 | |
1904 std::vector<const float*> rows(vertical.size()); | |
1905 | |
1906 for (unsigned int y = 0; y < height; y++) | |
1907 { | |
1908 for (unsigned int k = 0; k < vertical.size(); k++) | |
1909 { | |
1910 if (y < verticalAnchor) | |
1911 { | |
1912 rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(0)); // Use top border | |
1913 } | |
1914 else if (y + k >= height + verticalAnchor) | |
1915 { | |
1916 rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(height - 1)); // Use bottom border | |
1917 } | |
1918 else | |
1919 { | |
1920 rows[k] = reinterpret_cast<const float*>(tmp.GetConstRow(y + k - verticalAnchor)); | |
1921 } | |
1922 } | |
1923 | |
1924 RawPixel* p = reinterpret_cast<RawPixel*>(image.GetRow(y)); | |
1925 | |
1926 for (unsigned int x = 0; x < width; x++) | |
1927 { | |
1928 for (unsigned int c = 0; c < ChannelsCount; c++, p++) | |
1929 { | |
1930 float accumulator = 0; | |
1931 | |
1932 for (unsigned int k = 0; k < vertical.size(); k++) | |
1933 { | |
1934 accumulator += rows[k][ChannelsCount * x + c] * vertical[k]; | |
1935 } | |
1936 | |
1937 accumulator *= normalization; | |
1938 | |
1939 if (accumulator <= static_cast<float>(std::numeric_limits<RawPixel>::min())) | |
1940 { | |
1941 *p = std::numeric_limits<RawPixel>::min(); | |
1942 } | |
1943 else if (accumulator >= static_cast<float>(std::numeric_limits<RawPixel>::max())) | |
1944 { | |
1945 *p = std::numeric_limits<RawPixel>::max(); | |
1946 } | |
1947 else | |
1948 { | |
1949 *p = static_cast<RawPixel>(accumulator); | |
1950 } | |
1951 } | |
1952 } | |
1953 } | |
1954 } | |
1955 | |
1956 | |
1957 void ImageProcessing::SeparableConvolution(ImageAccessor& image /* inplace */, | |
1958 const std::vector<float>& horizontal, | |
1959 size_t horizontalAnchor, | |
1960 const std::vector<float>& vertical, | |
1961 size_t verticalAnchor) | |
1962 { | |
1963 if (horizontal.size() == 0 || | |
1964 vertical.size() == 0 || | |
1965 horizontalAnchor >= horizontal.size() || | |
1966 verticalAnchor >= vertical.size()) | |
1967 { | |
1968 throw OrthancException(ErrorCode_ParameterOutOfRange); | |
1969 } | |
1970 | |
1971 if (image.GetWidth() == 0 || | |
1972 image.GetHeight() == 0) | |
1973 { | |
1974 return; | |
1975 } | |
1976 | |
1977 /** | |
1978 * Compute normalization | |
1979 **/ | |
1980 | |
1981 float sumHorizontal = 0; | |
1982 for (size_t i = 0; i < horizontal.size(); i++) | |
1983 { | |
1984 sumHorizontal += horizontal[i]; | |
1985 } | |
1986 | |
1987 float sumVertical = 0; | |
1988 for (size_t i = 0; i < vertical.size(); i++) | |
1989 { | |
1990 sumVertical += vertical[i]; | |
1991 } | |
1992 | |
1993 if (fabsf(sumHorizontal) <= std::numeric_limits<float>::epsilon() || | |
1994 fabsf(sumVertical) <= std::numeric_limits<float>::epsilon()) | |
1995 { | |
1996 throw OrthancException(ErrorCode_ParameterOutOfRange, "Singular convolution kernel"); | |
1997 } | |
1998 | |
1999 const float normalization = 1.0f / (sumHorizontal * sumVertical); | |
2000 | |
2001 switch (image.GetFormat()) | |
2002 { | |
2003 case PixelFormat_Grayscale8: | |
2004 SeparableConvolutionFloat<uint8_t, 1u> | |
2005 (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization); | |
2006 break; | |
2007 | |
2008 case PixelFormat_RGB24: | |
2009 SeparableConvolutionFloat<uint8_t, 3u> | |
2010 (image, horizontal, horizontalAnchor, vertical, verticalAnchor, normalization); | |
2011 break; | |
2012 | |
2013 default: | |
2014 throw OrthancException(ErrorCode_NotImplemented); | |
2015 } | |
2016 } | |
2017 | |
2018 | |
2019 void ImageProcessing::SmoothGaussian5x5(ImageAccessor& image) | |
2020 { | |
2021 std::vector<float> kernel(5); | |
2022 kernel[0] = 1; | |
2023 kernel[1] = 4; | |
2024 kernel[2] = 6; | |
2025 kernel[3] = 4; | |
2026 kernel[4] = 1; | |
2027 | |
2028 SeparableConvolution(image, kernel, 2, kernel, 2); | |
2029 } | |
1768 } | 2030 } |