1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % FFFFF EEEEE AAA TTTTT U U RRRR EEEEE %
7 % F E A A T U U R R E %
8 % FFF EEE AAAAA T U U RRRR EEE %
9 % F E A A T U U R R E %
10 % F EEEEE A A T UUU R R EEEEE %
11 % %
12 % %
13 % MagickCore Image Feature Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright 1999-2021 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39
40 /*
41 Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/animate.h"
45 #include "MagickCore/artifact.h"
46 #include "MagickCore/blob.h"
47 #include "MagickCore/blob-private.h"
48 #include "MagickCore/cache.h"
49 #include "MagickCore/cache-private.h"
50 #include "MagickCore/cache-view.h"
51 #include "MagickCore/channel.h"
52 #include "MagickCore/client.h"
53 #include "MagickCore/color.h"
54 #include "MagickCore/color-private.h"
55 #include "MagickCore/colorspace.h"
56 #include "MagickCore/colorspace-private.h"
57 #include "MagickCore/composite.h"
58 #include "MagickCore/composite-private.h"
59 #include "MagickCore/compress.h"
60 #include "MagickCore/constitute.h"
61 #include "MagickCore/display.h"
62 #include "MagickCore/draw.h"
63 #include "MagickCore/enhance.h"
64 #include "MagickCore/exception.h"
65 #include "MagickCore/exception-private.h"
66 #include "MagickCore/feature.h"
67 #include "MagickCore/gem.h"
68 #include "MagickCore/geometry.h"
69 #include "MagickCore/list.h"
70 #include "MagickCore/image-private.h"
71 #include "MagickCore/magic.h"
72 #include "MagickCore/magick.h"
73 #include "MagickCore/matrix.h"
74 #include "MagickCore/memory_.h"
75 #include "MagickCore/module.h"
76 #include "MagickCore/monitor.h"
77 #include "MagickCore/monitor-private.h"
78 #include "MagickCore/morphology-private.h"
79 #include "MagickCore/option.h"
80 #include "MagickCore/paint.h"
81 #include "MagickCore/pixel-accessor.h"
82 #include "MagickCore/profile.h"
83 #include "MagickCore/property.h"
84 #include "MagickCore/quantize.h"
85 #include "MagickCore/quantum-private.h"
86 #include "MagickCore/random_.h"
87 #include "MagickCore/resource_.h"
88 #include "MagickCore/segment.h"
89 #include "MagickCore/semaphore.h"
90 #include "MagickCore/signature-private.h"
91 #include "MagickCore/string_.h"
92 #include "MagickCore/thread-private.h"
93 #include "MagickCore/timer.h"
94 #include "MagickCore/utility.h"
95 #include "MagickCore/version.h"
96
97 /*
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % %
100 % %
101 % %
102 % C a n n y E d g e I m a g e %
103 % %
104 % %
105 % %
106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 %
108 % CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
109 % edges in images.
110 %
111 % The format of the CannyEdgeImage method is:
112 %
113 % Image *CannyEdgeImage(const Image *image,const double radius,
114 % const double sigma,const double lower_percent,
115 % const double upper_percent,ExceptionInfo *exception)
116 %
117 % A description of each parameter follows:
118 %
119 % o image: the image.
120 %
121 % o radius: the radius of the gaussian smoothing filter.
122 %
123 % o sigma: the sigma of the gaussian smoothing filter.
124 %
125 % o lower_percent: percentage of edge pixels in the lower threshold.
126 %
127 % o upper_percent: percentage of edge pixels in the upper threshold.
128 %
129 % o exception: return any errors or warnings in this structure.
130 %
131 */
132
133 typedef struct _CannyInfo
134 {
135 double
136 magnitude,
137 intensity;
138
139 int
140 orientation;
141
142 ssize_t
143 x,
144 y;
145 } CannyInfo;
146
IsAuthenticPixel(const Image * image,const ssize_t x,const ssize_t y)147 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
148 const ssize_t x,const ssize_t y)
149 {
150 if ((x < 0) || (x >= (ssize_t) image->columns))
151 return(MagickFalse);
152 if ((y < 0) || (y >= (ssize_t) image->rows))
153 return(MagickFalse);
154 return(MagickTrue);
155 }
156
TraceEdges(Image * edge_image,CacheView * edge_view,MatrixInfo * canny_cache,const ssize_t x,const ssize_t y,const double lower_threshold,ExceptionInfo * exception)157 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
158 MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
159 const double lower_threshold,ExceptionInfo *exception)
160 {
161 CannyInfo
162 edge,
163 pixel;
164
165 MagickBooleanType
166 status;
167
168 Quantum
169 *q;
170
171 ssize_t
172 i;
173
174 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
175 if (q == (Quantum *) NULL)
176 return(MagickFalse);
177 *q=QuantumRange;
178 status=SyncCacheViewAuthenticPixels(edge_view,exception);
179 if (status == MagickFalse)
180 return(MagickFalse);
181 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
182 return(MagickFalse);
183 edge.x=x;
184 edge.y=y;
185 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
186 return(MagickFalse);
187 for (i=1; i != 0; )
188 {
189 ssize_t
190 v;
191
192 i--;
193 status=GetMatrixElement(canny_cache,i,0,&edge);
194 if (status == MagickFalse)
195 return(MagickFalse);
196 for (v=(-1); v <= 1; v++)
197 {
198 ssize_t
199 u;
200
201 for (u=(-1); u <= 1; u++)
202 {
203 if ((u == 0) && (v == 0))
204 continue;
205 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
206 continue;
207 /*
208 Not an edge if gradient value is below the lower threshold.
209 */
210 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
211 exception);
212 if (q == (Quantum *) NULL)
213 return(MagickFalse);
214 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
215 if (status == MagickFalse)
216 return(MagickFalse);
217 if ((GetPixelIntensity(edge_image,q) == 0.0) &&
218 (pixel.intensity >= lower_threshold))
219 {
220 *q=QuantumRange;
221 status=SyncCacheViewAuthenticPixels(edge_view,exception);
222 if (status == MagickFalse)
223 return(MagickFalse);
224 edge.x+=u;
225 edge.y+=v;
226 status=SetMatrixElement(canny_cache,i,0,&edge);
227 if (status == MagickFalse)
228 return(MagickFalse);
229 i++;
230 }
231 }
232 }
233 }
234 return(MagickTrue);
235 }
236
CannyEdgeImage(const Image * image,const double radius,const double sigma,const double lower_percent,const double upper_percent,ExceptionInfo * exception)237 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
238 const double sigma,const double lower_percent,const double upper_percent,
239 ExceptionInfo *exception)
240 {
241 #define CannyEdgeImageTag "CannyEdge/Image"
242
243 CacheView
244 *edge_view;
245
246 CannyInfo
247 element;
248
249 char
250 geometry[MagickPathExtent];
251
252 double
253 lower_threshold,
254 max,
255 min,
256 upper_threshold;
257
258 Image
259 *edge_image;
260
261 KernelInfo
262 *kernel_info;
263
264 MagickBooleanType
265 status;
266
267 MagickOffsetType
268 progress;
269
270 MatrixInfo
271 *canny_cache;
272
273 ssize_t
274 y;
275
276 assert(image != (const Image *) NULL);
277 assert(image->signature == MagickCoreSignature);
278 if (image->debug != MagickFalse)
279 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
280 assert(exception != (ExceptionInfo *) NULL);
281 assert(exception->signature == MagickCoreSignature);
282 /*
283 Filter out noise.
284 */
285 (void) FormatLocaleString(geometry,MagickPathExtent,
286 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
287 kernel_info=AcquireKernelInfo(geometry,exception);
288 if (kernel_info == (KernelInfo *) NULL)
289 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
290 edge_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,exception);
291 kernel_info=DestroyKernelInfo(kernel_info);
292 if (edge_image == (Image *) NULL)
293 return((Image *) NULL);
294 if (TransformImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
295 {
296 edge_image=DestroyImage(edge_image);
297 return((Image *) NULL);
298 }
299 (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
300 /*
301 Find the intensity gradient of the image.
302 */
303 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
304 sizeof(CannyInfo),exception);
305 if (canny_cache == (MatrixInfo *) NULL)
306 {
307 edge_image=DestroyImage(edge_image);
308 return((Image *) NULL);
309 }
310 status=MagickTrue;
311 edge_view=AcquireVirtualCacheView(edge_image,exception);
312 #if defined(MAGICKCORE_OPENMP_SUPPORT)
313 #pragma omp parallel for schedule(static) shared(status) \
314 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
315 #endif
316 for (y=0; y < (ssize_t) edge_image->rows; y++)
317 {
318 const Quantum
319 *magick_restrict p;
320
321 ssize_t
322 x;
323
324 if (status == MagickFalse)
325 continue;
326 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
327 exception);
328 if (p == (const Quantum *) NULL)
329 {
330 status=MagickFalse;
331 continue;
332 }
333 for (x=0; x < (ssize_t) edge_image->columns; x++)
334 {
335 CannyInfo
336 pixel;
337
338 double
339 dx,
340 dy;
341
342 const Quantum
343 *magick_restrict kernel_pixels;
344
345 ssize_t
346 v;
347
348 static double
349 Gx[2][2] =
350 {
351 { -1.0, +1.0 },
352 { -1.0, +1.0 }
353 },
354 Gy[2][2] =
355 {
356 { +1.0, +1.0 },
357 { -1.0, -1.0 }
358 };
359
360 (void) memset(&pixel,0,sizeof(pixel));
361 dx=0.0;
362 dy=0.0;
363 kernel_pixels=p;
364 for (v=0; v < 2; v++)
365 {
366 ssize_t
367 u;
368
369 for (u=0; u < 2; u++)
370 {
371 double
372 intensity;
373
374 intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
375 dx+=0.5*Gx[v][u]*intensity;
376 dy+=0.5*Gy[v][u]*intensity;
377 }
378 kernel_pixels+=edge_image->columns+1;
379 }
380 pixel.magnitude=hypot(dx,dy);
381 pixel.orientation=0;
382 if (fabs(dx) > MagickEpsilon)
383 {
384 double
385 slope;
386
387 slope=dy/dx;
388 if (slope < 0.0)
389 {
390 if (slope < -2.41421356237)
391 pixel.orientation=0;
392 else
393 if (slope < -0.414213562373)
394 pixel.orientation=1;
395 else
396 pixel.orientation=2;
397 }
398 else
399 {
400 if (slope > 2.41421356237)
401 pixel.orientation=0;
402 else
403 if (slope > 0.414213562373)
404 pixel.orientation=3;
405 else
406 pixel.orientation=2;
407 }
408 }
409 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
410 continue;
411 p+=GetPixelChannels(edge_image);
412 }
413 }
414 edge_view=DestroyCacheView(edge_view);
415 /*
416 Non-maxima suppression, remove pixels that are not considered to be part
417 of an edge.
418 */
419 progress=0;
420 (void) GetMatrixElement(canny_cache,0,0,&element);
421 max=element.intensity;
422 min=element.intensity;
423 edge_view=AcquireAuthenticCacheView(edge_image,exception);
424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
425 #pragma omp parallel for schedule(static) shared(status) \
426 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
427 #endif
428 for (y=0; y < (ssize_t) edge_image->rows; y++)
429 {
430 Quantum
431 *magick_restrict q;
432
433 ssize_t
434 x;
435
436 if (status == MagickFalse)
437 continue;
438 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
439 exception);
440 if (q == (Quantum *) NULL)
441 {
442 status=MagickFalse;
443 continue;
444 }
445 for (x=0; x < (ssize_t) edge_image->columns; x++)
446 {
447 CannyInfo
448 alpha_pixel,
449 beta_pixel,
450 pixel;
451
452 (void) GetMatrixElement(canny_cache,x,y,&pixel);
453 switch (pixel.orientation)
454 {
455 case 0:
456 default:
457 {
458 /*
459 0 degrees, north and south.
460 */
461 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
462 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
463 break;
464 }
465 case 1:
466 {
467 /*
468 45 degrees, northwest and southeast.
469 */
470 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
471 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
472 break;
473 }
474 case 2:
475 {
476 /*
477 90 degrees, east and west.
478 */
479 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
480 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
481 break;
482 }
483 case 3:
484 {
485 /*
486 135 degrees, northeast and southwest.
487 */
488 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
489 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
490 break;
491 }
492 }
493 pixel.intensity=pixel.magnitude;
494 if ((pixel.magnitude < alpha_pixel.magnitude) ||
495 (pixel.magnitude < beta_pixel.magnitude))
496 pixel.intensity=0;
497 (void) SetMatrixElement(canny_cache,x,y,&pixel);
498 #if defined(MAGICKCORE_OPENMP_SUPPORT)
499 #pragma omp critical (MagickCore_CannyEdgeImage)
500 #endif
501 {
502 if (pixel.intensity < min)
503 min=pixel.intensity;
504 if (pixel.intensity > max)
505 max=pixel.intensity;
506 }
507 *q=0;
508 q+=GetPixelChannels(edge_image);
509 }
510 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
511 status=MagickFalse;
512 }
513 edge_view=DestroyCacheView(edge_view);
514 /*
515 Estimate hysteresis threshold.
516 */
517 lower_threshold=lower_percent*(max-min)+min;
518 upper_threshold=upper_percent*(max-min)+min;
519 /*
520 Hysteresis threshold.
521 */
522 edge_view=AcquireAuthenticCacheView(edge_image,exception);
523 for (y=0; y < (ssize_t) edge_image->rows; y++)
524 {
525 ssize_t
526 x;
527
528 if (status == MagickFalse)
529 continue;
530 for (x=0; x < (ssize_t) edge_image->columns; x++)
531 {
532 CannyInfo
533 pixel;
534
535 const Quantum
536 *magick_restrict p;
537
538 /*
539 Edge if pixel gradient higher than upper threshold.
540 */
541 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
542 if (p == (const Quantum *) NULL)
543 continue;
544 status=GetMatrixElement(canny_cache,x,y,&pixel);
545 if (status == MagickFalse)
546 continue;
547 if ((GetPixelIntensity(edge_image,p) == 0.0) &&
548 (pixel.intensity >= upper_threshold))
549 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
550 exception);
551 }
552 if (image->progress_monitor != (MagickProgressMonitor) NULL)
553 {
554 MagickBooleanType
555 proceed;
556
557 #if defined(MAGICKCORE_OPENMP_SUPPORT)
558 #pragma omp atomic
559 #endif
560 progress++;
561 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
562 if (proceed == MagickFalse)
563 status=MagickFalse;
564 }
565 }
566 edge_view=DestroyCacheView(edge_view);
567 /*
568 Free resources.
569 */
570 canny_cache=DestroyMatrixInfo(canny_cache);
571 return(edge_image);
572 }
573
574 /*
575 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
576 % %
577 % %
578 % %
579 % G e t I m a g e F e a t u r e s %
580 % %
581 % %
582 % %
583 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
584 %
585 % GetImageFeatures() returns features for each channel in the image in
586 % each of four directions (horizontal, vertical, left and right diagonals)
587 % for the specified distance. The features include the angular second
588 % moment, contrast, correlation, sum of squares: variance, inverse difference
589 % moment, sum average, sum varience, sum entropy, entropy, difference variance,
590 % difference entropy, information measures of correlation 1, information
591 % measures of correlation 2, and maximum correlation coefficient. You can
592 % access the red channel contrast, for example, like this:
593 %
594 % channel_features=GetImageFeatures(image,1,exception);
595 % contrast=channel_features[RedPixelChannel].contrast[0];
596 %
597 % Use MagickRelinquishMemory() to free the features buffer.
598 %
599 % The format of the GetImageFeatures method is:
600 %
601 % ChannelFeatures *GetImageFeatures(const Image *image,
602 % const size_t distance,ExceptionInfo *exception)
603 %
604 % A description of each parameter follows:
605 %
606 % o image: the image.
607 %
608 % o distance: the distance.
609 %
610 % o exception: return any errors or warnings in this structure.
611 %
612 */
613
MagickLog10(const double x)614 static inline double MagickLog10(const double x)
615 {
616 #define Log10Epsilon (1.0e-11)
617
618 if (fabs(x) < Log10Epsilon)
619 return(log10(Log10Epsilon));
620 return(log10(fabs(x)));
621 }
622
GetImageFeatures(const Image * image,const size_t distance,ExceptionInfo * exception)623 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
624 const size_t distance,ExceptionInfo *exception)
625 {
626 typedef struct _ChannelStatistics
627 {
628 PixelInfo
629 direction[4]; /* horizontal, vertical, left and right diagonals */
630 } ChannelStatistics;
631
632 CacheView
633 *image_view;
634
635 ChannelFeatures
636 *channel_features;
637
638 ChannelStatistics
639 **cooccurrence,
640 correlation,
641 *density_x,
642 *density_xy,
643 *density_y,
644 entropy_x,
645 entropy_xy,
646 entropy_xy1,
647 entropy_xy2,
648 entropy_y,
649 mean,
650 **Q,
651 *sum,
652 sum_squares,
653 variance;
654
655 PixelPacket
656 gray,
657 *grays;
658
659 MagickBooleanType
660 status;
661
662 ssize_t
663 i,
664 r;
665
666 size_t
667 length;
668
669 unsigned int
670 number_grays;
671
672 assert(image != (Image *) NULL);
673 assert(image->signature == MagickCoreSignature);
674 if (image->debug != MagickFalse)
675 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
676 if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
677 return((ChannelFeatures *) NULL);
678 length=MaxPixelChannels+1UL;
679 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
680 sizeof(*channel_features));
681 if (channel_features == (ChannelFeatures *) NULL)
682 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
683 (void) memset(channel_features,0,length*
684 sizeof(*channel_features));
685 /*
686 Form grays.
687 */
688 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
689 if (grays == (PixelPacket *) NULL)
690 {
691 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
692 channel_features);
693 (void) ThrowMagickException(exception,GetMagickModule(),
694 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
695 return(channel_features);
696 }
697 for (i=0; i <= (ssize_t) MaxMap; i++)
698 {
699 grays[i].red=(~0U);
700 grays[i].green=(~0U);
701 grays[i].blue=(~0U);
702 grays[i].alpha=(~0U);
703 grays[i].black=(~0U);
704 }
705 status=MagickTrue;
706 image_view=AcquireVirtualCacheView(image,exception);
707 #if defined(MAGICKCORE_OPENMP_SUPPORT)
708 #pragma omp parallel for schedule(static) shared(status) \
709 magick_number_threads(image,image,image->rows,1)
710 #endif
711 for (r=0; r < (ssize_t) image->rows; r++)
712 {
713 const Quantum
714 *magick_restrict p;
715
716 ssize_t
717 x;
718
719 if (status == MagickFalse)
720 continue;
721 p=GetCacheViewVirtualPixels(image_view,0,r,image->columns,1,exception);
722 if (p == (const Quantum *) NULL)
723 {
724 status=MagickFalse;
725 continue;
726 }
727 for (x=0; x < (ssize_t) image->columns; x++)
728 {
729 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
730 ScaleQuantumToMap(GetPixelRed(image,p));
731 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
732 ScaleQuantumToMap(GetPixelGreen(image,p));
733 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
734 ScaleQuantumToMap(GetPixelBlue(image,p));
735 if (image->colorspace == CMYKColorspace)
736 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
737 ScaleQuantumToMap(GetPixelBlack(image,p));
738 if (image->alpha_trait != UndefinedPixelTrait)
739 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
740 ScaleQuantumToMap(GetPixelAlpha(image,p));
741 p+=GetPixelChannels(image);
742 }
743 }
744 image_view=DestroyCacheView(image_view);
745 if (status == MagickFalse)
746 {
747 grays=(PixelPacket *) RelinquishMagickMemory(grays);
748 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
749 channel_features);
750 return(channel_features);
751 }
752 (void) memset(&gray,0,sizeof(gray));
753 for (i=0; i <= (ssize_t) MaxMap; i++)
754 {
755 if (grays[i].red != ~0U)
756 grays[gray.red++].red=grays[i].red;
757 if (grays[i].green != ~0U)
758 grays[gray.green++].green=grays[i].green;
759 if (grays[i].blue != ~0U)
760 grays[gray.blue++].blue=grays[i].blue;
761 if (image->colorspace == CMYKColorspace)
762 if (grays[i].black != ~0U)
763 grays[gray.black++].black=grays[i].black;
764 if (image->alpha_trait != UndefinedPixelTrait)
765 if (grays[i].alpha != ~0U)
766 grays[gray.alpha++].alpha=grays[i].alpha;
767 }
768 /*
769 Allocate spatial dependence matrix.
770 */
771 number_grays=gray.red;
772 if (gray.green > number_grays)
773 number_grays=gray.green;
774 if (gray.blue > number_grays)
775 number_grays=gray.blue;
776 if (image->colorspace == CMYKColorspace)
777 if (gray.black > number_grays)
778 number_grays=gray.black;
779 if (image->alpha_trait != UndefinedPixelTrait)
780 if (gray.alpha > number_grays)
781 number_grays=gray.alpha;
782 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
783 sizeof(*cooccurrence));
784 density_x=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
785 2*sizeof(*density_x));
786 density_xy=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
787 2*sizeof(*density_xy));
788 density_y=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
789 2*sizeof(*density_y));
790 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
791 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
792 if ((cooccurrence == (ChannelStatistics **) NULL) ||
793 (density_x == (ChannelStatistics *) NULL) ||
794 (density_xy == (ChannelStatistics *) NULL) ||
795 (density_y == (ChannelStatistics *) NULL) ||
796 (Q == (ChannelStatistics **) NULL) ||
797 (sum == (ChannelStatistics *) NULL))
798 {
799 if (Q != (ChannelStatistics **) NULL)
800 {
801 for (i=0; i < (ssize_t) number_grays; i++)
802 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
803 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
804 }
805 if (sum != (ChannelStatistics *) NULL)
806 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
807 if (density_y != (ChannelStatistics *) NULL)
808 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
809 if (density_xy != (ChannelStatistics *) NULL)
810 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
811 if (density_x != (ChannelStatistics *) NULL)
812 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
813 if (cooccurrence != (ChannelStatistics **) NULL)
814 {
815 for (i=0; i < (ssize_t) number_grays; i++)
816 cooccurrence[i]=(ChannelStatistics *)
817 RelinquishMagickMemory(cooccurrence[i]);
818 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
819 cooccurrence);
820 }
821 grays=(PixelPacket *) RelinquishMagickMemory(grays);
822 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
823 channel_features);
824 (void) ThrowMagickException(exception,GetMagickModule(),
825 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
826 return(channel_features);
827 }
828 (void) memset(&correlation,0,sizeof(correlation));
829 (void) memset(density_x,0,2*(number_grays+1)*sizeof(*density_x));
830 (void) memset(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
831 (void) memset(density_y,0,2*(number_grays+1)*sizeof(*density_y));
832 (void) memset(&mean,0,sizeof(mean));
833 (void) memset(sum,0,number_grays*sizeof(*sum));
834 (void) memset(&sum_squares,0,sizeof(sum_squares));
835 (void) memset(density_xy,0,2*number_grays*sizeof(*density_xy));
836 (void) memset(&entropy_x,0,sizeof(entropy_x));
837 (void) memset(&entropy_xy,0,sizeof(entropy_xy));
838 (void) memset(&entropy_xy1,0,sizeof(entropy_xy1));
839 (void) memset(&entropy_xy2,0,sizeof(entropy_xy2));
840 (void) memset(&entropy_y,0,sizeof(entropy_y));
841 (void) memset(&variance,0,sizeof(variance));
842 for (i=0; i < (ssize_t) number_grays; i++)
843 {
844 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
845 sizeof(**cooccurrence));
846 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
847 if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
848 (Q[i] == (ChannelStatistics *) NULL))
849 break;
850 (void) memset(cooccurrence[i],0,number_grays*
851 sizeof(**cooccurrence));
852 (void) memset(Q[i],0,number_grays*sizeof(**Q));
853 }
854 if (i < (ssize_t) number_grays)
855 {
856 for (i--; i >= 0; i--)
857 {
858 if (Q[i] != (ChannelStatistics *) NULL)
859 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
860 if (cooccurrence[i] != (ChannelStatistics *) NULL)
861 cooccurrence[i]=(ChannelStatistics *)
862 RelinquishMagickMemory(cooccurrence[i]);
863 }
864 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
865 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
866 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
867 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
868 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
869 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
870 grays=(PixelPacket *) RelinquishMagickMemory(grays);
871 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
872 channel_features);
873 (void) ThrowMagickException(exception,GetMagickModule(),
874 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
875 return(channel_features);
876 }
877 /*
878 Initialize spatial dependence matrix.
879 */
880 status=MagickTrue;
881 image_view=AcquireVirtualCacheView(image,exception);
882 for (r=0; r < (ssize_t) image->rows; r++)
883 {
884 const Quantum
885 *magick_restrict p;
886
887 ssize_t
888 x;
889
890 ssize_t
891 offset,
892 u,
893 v;
894
895 if (status == MagickFalse)
896 continue;
897 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+
898 2*distance,distance+2,exception);
899 if (p == (const Quantum *) NULL)
900 {
901 status=MagickFalse;
902 continue;
903 }
904 p+=distance*GetPixelChannels(image);;
905 for (x=0; x < (ssize_t) image->columns; x++)
906 {
907 for (i=0; i < 4; i++)
908 {
909 switch (i)
910 {
911 case 0:
912 default:
913 {
914 /*
915 Horizontal adjacency.
916 */
917 offset=(ssize_t) distance;
918 break;
919 }
920 case 1:
921 {
922 /*
923 Vertical adjacency.
924 */
925 offset=(ssize_t) (image->columns+2*distance);
926 break;
927 }
928 case 2:
929 {
930 /*
931 Right diagonal adjacency.
932 */
933 offset=(ssize_t) ((image->columns+2*distance)-distance);
934 break;
935 }
936 case 3:
937 {
938 /*
939 Left diagonal adjacency.
940 */
941 offset=(ssize_t) ((image->columns+2*distance)+distance);
942 break;
943 }
944 }
945 u=0;
946 v=0;
947 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
948 u++;
949 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
950 v++;
951 cooccurrence[u][v].direction[i].red++;
952 cooccurrence[v][u].direction[i].red++;
953 u=0;
954 v=0;
955 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
956 u++;
957 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
958 v++;
959 cooccurrence[u][v].direction[i].green++;
960 cooccurrence[v][u].direction[i].green++;
961 u=0;
962 v=0;
963 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
964 u++;
965 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
966 v++;
967 cooccurrence[u][v].direction[i].blue++;
968 cooccurrence[v][u].direction[i].blue++;
969 if (image->colorspace == CMYKColorspace)
970 {
971 u=0;
972 v=0;
973 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
974 u++;
975 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
976 v++;
977 cooccurrence[u][v].direction[i].black++;
978 cooccurrence[v][u].direction[i].black++;
979 }
980 if (image->alpha_trait != UndefinedPixelTrait)
981 {
982 u=0;
983 v=0;
984 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
985 u++;
986 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
987 v++;
988 cooccurrence[u][v].direction[i].alpha++;
989 cooccurrence[v][u].direction[i].alpha++;
990 }
991 }
992 p+=GetPixelChannels(image);
993 }
994 }
995 grays=(PixelPacket *) RelinquishMagickMemory(grays);
996 image_view=DestroyCacheView(image_view);
997 if (status == MagickFalse)
998 {
999 for (i=0; i < (ssize_t) number_grays; i++)
1000 cooccurrence[i]=(ChannelStatistics *)
1001 RelinquishMagickMemory(cooccurrence[i]);
1002 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1003 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1004 channel_features);
1005 (void) ThrowMagickException(exception,GetMagickModule(),
1006 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1007 return(channel_features);
1008 }
1009 /*
1010 Normalize spatial dependence matrix.
1011 */
1012 for (i=0; i < 4; i++)
1013 {
1014 double
1015 normalize;
1016
1017 ssize_t
1018 y;
1019
1020 switch (i)
1021 {
1022 case 0:
1023 default:
1024 {
1025 /*
1026 Horizontal adjacency.
1027 */
1028 normalize=2.0*image->rows*(image->columns-distance);
1029 break;
1030 }
1031 case 1:
1032 {
1033 /*
1034 Vertical adjacency.
1035 */
1036 normalize=2.0*(image->rows-distance)*image->columns;
1037 break;
1038 }
1039 case 2:
1040 {
1041 /*
1042 Right diagonal adjacency.
1043 */
1044 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1045 break;
1046 }
1047 case 3:
1048 {
1049 /*
1050 Left diagonal adjacency.
1051 */
1052 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1053 break;
1054 }
1055 }
1056 normalize=PerceptibleReciprocal(normalize);
1057 for (y=0; y < (ssize_t) number_grays; y++)
1058 {
1059 ssize_t
1060 x;
1061
1062 for (x=0; x < (ssize_t) number_grays; x++)
1063 {
1064 cooccurrence[x][y].direction[i].red*=normalize;
1065 cooccurrence[x][y].direction[i].green*=normalize;
1066 cooccurrence[x][y].direction[i].blue*=normalize;
1067 if (image->colorspace == CMYKColorspace)
1068 cooccurrence[x][y].direction[i].black*=normalize;
1069 if (image->alpha_trait != UndefinedPixelTrait)
1070 cooccurrence[x][y].direction[i].alpha*=normalize;
1071 }
1072 }
1073 }
1074 /*
1075 Compute texture features.
1076 */
1077 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1078 #pragma omp parallel for schedule(static) shared(status) \
1079 magick_number_threads(image,image,number_grays,1)
1080 #endif
1081 for (i=0; i < 4; i++)
1082 {
1083 ssize_t
1084 y;
1085
1086 for (y=0; y < (ssize_t) number_grays; y++)
1087 {
1088 ssize_t
1089 x;
1090
1091 for (x=0; x < (ssize_t) number_grays; x++)
1092 {
1093 /*
1094 Angular second moment: measure of homogeneity of the image.
1095 */
1096 channel_features[RedPixelChannel].angular_second_moment[i]+=
1097 cooccurrence[x][y].direction[i].red*
1098 cooccurrence[x][y].direction[i].red;
1099 channel_features[GreenPixelChannel].angular_second_moment[i]+=
1100 cooccurrence[x][y].direction[i].green*
1101 cooccurrence[x][y].direction[i].green;
1102 channel_features[BluePixelChannel].angular_second_moment[i]+=
1103 cooccurrence[x][y].direction[i].blue*
1104 cooccurrence[x][y].direction[i].blue;
1105 if (image->colorspace == CMYKColorspace)
1106 channel_features[BlackPixelChannel].angular_second_moment[i]+=
1107 cooccurrence[x][y].direction[i].black*
1108 cooccurrence[x][y].direction[i].black;
1109 if (image->alpha_trait != UndefinedPixelTrait)
1110 channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1111 cooccurrence[x][y].direction[i].alpha*
1112 cooccurrence[x][y].direction[i].alpha;
1113 /*
1114 Correlation: measure of linear-dependencies in the image.
1115 */
1116 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1117 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1118 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1119 if (image->colorspace == CMYKColorspace)
1120 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1121 if (image->alpha_trait != UndefinedPixelTrait)
1122 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1123 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1124 correlation.direction[i].green+=x*y*
1125 cooccurrence[x][y].direction[i].green;
1126 correlation.direction[i].blue+=x*y*
1127 cooccurrence[x][y].direction[i].blue;
1128 if (image->colorspace == CMYKColorspace)
1129 correlation.direction[i].black+=x*y*
1130 cooccurrence[x][y].direction[i].black;
1131 if (image->alpha_trait != UndefinedPixelTrait)
1132 correlation.direction[i].alpha+=x*y*
1133 cooccurrence[x][y].direction[i].alpha;
1134 /*
1135 Inverse Difference Moment.
1136 */
1137 channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1138 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1139 channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1140 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1141 channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1142 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1143 if (image->colorspace == CMYKColorspace)
1144 channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1145 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1146 if (image->alpha_trait != UndefinedPixelTrait)
1147 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1148 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1149 /*
1150 Sum average.
1151 */
1152 density_xy[y+x+2].direction[i].red+=
1153 cooccurrence[x][y].direction[i].red;
1154 density_xy[y+x+2].direction[i].green+=
1155 cooccurrence[x][y].direction[i].green;
1156 density_xy[y+x+2].direction[i].blue+=
1157 cooccurrence[x][y].direction[i].blue;
1158 if (image->colorspace == CMYKColorspace)
1159 density_xy[y+x+2].direction[i].black+=
1160 cooccurrence[x][y].direction[i].black;
1161 if (image->alpha_trait != UndefinedPixelTrait)
1162 density_xy[y+x+2].direction[i].alpha+=
1163 cooccurrence[x][y].direction[i].alpha;
1164 /*
1165 Entropy.
1166 */
1167 channel_features[RedPixelChannel].entropy[i]-=
1168 cooccurrence[x][y].direction[i].red*
1169 MagickLog10(cooccurrence[x][y].direction[i].red);
1170 channel_features[GreenPixelChannel].entropy[i]-=
1171 cooccurrence[x][y].direction[i].green*
1172 MagickLog10(cooccurrence[x][y].direction[i].green);
1173 channel_features[BluePixelChannel].entropy[i]-=
1174 cooccurrence[x][y].direction[i].blue*
1175 MagickLog10(cooccurrence[x][y].direction[i].blue);
1176 if (image->colorspace == CMYKColorspace)
1177 channel_features[BlackPixelChannel].entropy[i]-=
1178 cooccurrence[x][y].direction[i].black*
1179 MagickLog10(cooccurrence[x][y].direction[i].black);
1180 if (image->alpha_trait != UndefinedPixelTrait)
1181 channel_features[AlphaPixelChannel].entropy[i]-=
1182 cooccurrence[x][y].direction[i].alpha*
1183 MagickLog10(cooccurrence[x][y].direction[i].alpha);
1184 /*
1185 Information Measures of Correlation.
1186 */
1187 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1188 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1189 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1190 if (image->alpha_trait != UndefinedPixelTrait)
1191 density_x[x].direction[i].alpha+=
1192 cooccurrence[x][y].direction[i].alpha;
1193 if (image->colorspace == CMYKColorspace)
1194 density_x[x].direction[i].black+=
1195 cooccurrence[x][y].direction[i].black;
1196 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1197 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1198 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1199 if (image->colorspace == CMYKColorspace)
1200 density_y[y].direction[i].black+=
1201 cooccurrence[x][y].direction[i].black;
1202 if (image->alpha_trait != UndefinedPixelTrait)
1203 density_y[y].direction[i].alpha+=
1204 cooccurrence[x][y].direction[i].alpha;
1205 }
1206 mean.direction[i].red+=y*sum[y].direction[i].red;
1207 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1208 mean.direction[i].green+=y*sum[y].direction[i].green;
1209 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1210 mean.direction[i].blue+=y*sum[y].direction[i].blue;
1211 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1212 if (image->colorspace == CMYKColorspace)
1213 {
1214 mean.direction[i].black+=y*sum[y].direction[i].black;
1215 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1216 }
1217 if (image->alpha_trait != UndefinedPixelTrait)
1218 {
1219 mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1220 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1221 }
1222 }
1223 /*
1224 Correlation: measure of linear-dependencies in the image.
1225 */
1226 channel_features[RedPixelChannel].correlation[i]=
1227 (correlation.direction[i].red-mean.direction[i].red*
1228 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1229 (mean.direction[i].red*mean.direction[i].red))*sqrt(
1230 sum_squares.direction[i].red-(mean.direction[i].red*
1231 mean.direction[i].red)));
1232 channel_features[GreenPixelChannel].correlation[i]=
1233 (correlation.direction[i].green-mean.direction[i].green*
1234 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1235 (mean.direction[i].green*mean.direction[i].green))*sqrt(
1236 sum_squares.direction[i].green-(mean.direction[i].green*
1237 mean.direction[i].green)));
1238 channel_features[BluePixelChannel].correlation[i]=
1239 (correlation.direction[i].blue-mean.direction[i].blue*
1240 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1241 (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1242 sum_squares.direction[i].blue-(mean.direction[i].blue*
1243 mean.direction[i].blue)));
1244 if (image->colorspace == CMYKColorspace)
1245 channel_features[BlackPixelChannel].correlation[i]=
1246 (correlation.direction[i].black-mean.direction[i].black*
1247 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1248 (mean.direction[i].black*mean.direction[i].black))*sqrt(
1249 sum_squares.direction[i].black-(mean.direction[i].black*
1250 mean.direction[i].black)));
1251 if (image->alpha_trait != UndefinedPixelTrait)
1252 channel_features[AlphaPixelChannel].correlation[i]=
1253 (correlation.direction[i].alpha-mean.direction[i].alpha*
1254 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1255 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1256 sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1257 mean.direction[i].alpha)));
1258 }
1259 /*
1260 Compute more texture features.
1261 */
1262 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1263 #pragma omp parallel for schedule(static) shared(status) \
1264 magick_number_threads(image,image,number_grays,1)
1265 #endif
1266 for (i=0; i < 4; i++)
1267 {
1268 ssize_t
1269 x;
1270
1271 for (x=2; x < (ssize_t) (2*number_grays); x++)
1272 {
1273 /*
1274 Sum average.
1275 */
1276 channel_features[RedPixelChannel].sum_average[i]+=
1277 x*density_xy[x].direction[i].red;
1278 channel_features[GreenPixelChannel].sum_average[i]+=
1279 x*density_xy[x].direction[i].green;
1280 channel_features[BluePixelChannel].sum_average[i]+=
1281 x*density_xy[x].direction[i].blue;
1282 if (image->colorspace == CMYKColorspace)
1283 channel_features[BlackPixelChannel].sum_average[i]+=
1284 x*density_xy[x].direction[i].black;
1285 if (image->alpha_trait != UndefinedPixelTrait)
1286 channel_features[AlphaPixelChannel].sum_average[i]+=
1287 x*density_xy[x].direction[i].alpha;
1288 /*
1289 Sum entropy.
1290 */
1291 channel_features[RedPixelChannel].sum_entropy[i]-=
1292 density_xy[x].direction[i].red*
1293 MagickLog10(density_xy[x].direction[i].red);
1294 channel_features[GreenPixelChannel].sum_entropy[i]-=
1295 density_xy[x].direction[i].green*
1296 MagickLog10(density_xy[x].direction[i].green);
1297 channel_features[BluePixelChannel].sum_entropy[i]-=
1298 density_xy[x].direction[i].blue*
1299 MagickLog10(density_xy[x].direction[i].blue);
1300 if (image->colorspace == CMYKColorspace)
1301 channel_features[BlackPixelChannel].sum_entropy[i]-=
1302 density_xy[x].direction[i].black*
1303 MagickLog10(density_xy[x].direction[i].black);
1304 if (image->alpha_trait != UndefinedPixelTrait)
1305 channel_features[AlphaPixelChannel].sum_entropy[i]-=
1306 density_xy[x].direction[i].alpha*
1307 MagickLog10(density_xy[x].direction[i].alpha);
1308 /*
1309 Sum variance.
1310 */
1311 channel_features[RedPixelChannel].sum_variance[i]+=
1312 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1313 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1314 density_xy[x].direction[i].red;
1315 channel_features[GreenPixelChannel].sum_variance[i]+=
1316 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1317 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1318 density_xy[x].direction[i].green;
1319 channel_features[BluePixelChannel].sum_variance[i]+=
1320 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1321 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1322 density_xy[x].direction[i].blue;
1323 if (image->colorspace == CMYKColorspace)
1324 channel_features[BlackPixelChannel].sum_variance[i]+=
1325 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1326 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1327 density_xy[x].direction[i].black;
1328 if (image->alpha_trait != UndefinedPixelTrait)
1329 channel_features[AlphaPixelChannel].sum_variance[i]+=
1330 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1331 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1332 density_xy[x].direction[i].alpha;
1333 }
1334 }
1335 /*
1336 Compute more texture features.
1337 */
1338 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1339 #pragma omp parallel for schedule(static) shared(status) \
1340 magick_number_threads(image,image,number_grays,1)
1341 #endif
1342 for (i=0; i < 4; i++)
1343 {
1344 ssize_t
1345 y;
1346
1347 for (y=0; y < (ssize_t) number_grays; y++)
1348 {
1349 ssize_t
1350 x;
1351
1352 for (x=0; x < (ssize_t) number_grays; x++)
1353 {
1354 /*
1355 Sum of Squares: Variance
1356 */
1357 variance.direction[i].red+=(y-mean.direction[i].red+1)*
1358 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1359 variance.direction[i].green+=(y-mean.direction[i].green+1)*
1360 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1361 variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1362 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1363 if (image->colorspace == CMYKColorspace)
1364 variance.direction[i].black+=(y-mean.direction[i].black+1)*
1365 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1366 if (image->alpha_trait != UndefinedPixelTrait)
1367 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1368 (y-mean.direction[i].alpha+1)*
1369 cooccurrence[x][y].direction[i].alpha;
1370 /*
1371 Sum average / Difference Variance.
1372 */
1373 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1374 cooccurrence[x][y].direction[i].red;
1375 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1376 cooccurrence[x][y].direction[i].green;
1377 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1378 cooccurrence[x][y].direction[i].blue;
1379 if (image->colorspace == CMYKColorspace)
1380 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1381 cooccurrence[x][y].direction[i].black;
1382 if (image->alpha_trait != UndefinedPixelTrait)
1383 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1384 cooccurrence[x][y].direction[i].alpha;
1385 /*
1386 Information Measures of Correlation.
1387 */
1388 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1389 MagickLog10(cooccurrence[x][y].direction[i].red);
1390 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1391 MagickLog10(cooccurrence[x][y].direction[i].green);
1392 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1393 MagickLog10(cooccurrence[x][y].direction[i].blue);
1394 if (image->colorspace == CMYKColorspace)
1395 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1396 MagickLog10(cooccurrence[x][y].direction[i].black);
1397 if (image->alpha_trait != UndefinedPixelTrait)
1398 entropy_xy.direction[i].alpha-=
1399 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1400 cooccurrence[x][y].direction[i].alpha);
1401 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1402 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1403 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1404 MagickLog10(density_x[x].direction[i].green*
1405 density_y[y].direction[i].green));
1406 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1407 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1408 if (image->colorspace == CMYKColorspace)
1409 entropy_xy1.direction[i].black-=(
1410 cooccurrence[x][y].direction[i].black*MagickLog10(
1411 density_x[x].direction[i].black*density_y[y].direction[i].black));
1412 if (image->alpha_trait != UndefinedPixelTrait)
1413 entropy_xy1.direction[i].alpha-=(
1414 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1415 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1416 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1417 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1418 density_y[y].direction[i].red));
1419 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1420 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1421 density_y[y].direction[i].green));
1422 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1423 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1424 density_y[y].direction[i].blue));
1425 if (image->colorspace == CMYKColorspace)
1426 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1427 density_y[y].direction[i].black*MagickLog10(
1428 density_x[x].direction[i].black*density_y[y].direction[i].black));
1429 if (image->alpha_trait != UndefinedPixelTrait)
1430 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1431 density_y[y].direction[i].alpha*MagickLog10(
1432 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1433 }
1434 }
1435 channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1436 variance.direction[i].red;
1437 channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1438 variance.direction[i].green;
1439 channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1440 variance.direction[i].blue;
1441 if (image->colorspace == CMYKColorspace)
1442 channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1443 variance.direction[i].black;
1444 if (image->alpha_trait != UndefinedPixelTrait)
1445 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1446 variance.direction[i].alpha;
1447 }
1448 /*
1449 Compute more texture features.
1450 */
1451 (void) memset(&variance,0,sizeof(variance));
1452 (void) memset(&sum_squares,0,sizeof(sum_squares));
1453 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1454 #pragma omp parallel for schedule(static) shared(status) \
1455 magick_number_threads(image,image,number_grays,1)
1456 #endif
1457 for (i=0; i < 4; i++)
1458 {
1459 ssize_t
1460 x;
1461
1462 for (x=0; x < (ssize_t) number_grays; x++)
1463 {
1464 /*
1465 Difference variance.
1466 */
1467 variance.direction[i].red+=density_xy[x].direction[i].red;
1468 variance.direction[i].green+=density_xy[x].direction[i].green;
1469 variance.direction[i].blue+=density_xy[x].direction[i].blue;
1470 if (image->colorspace == CMYKColorspace)
1471 variance.direction[i].black+=density_xy[x].direction[i].black;
1472 if (image->alpha_trait != UndefinedPixelTrait)
1473 variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1474 sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1475 density_xy[x].direction[i].red;
1476 sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1477 density_xy[x].direction[i].green;
1478 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1479 density_xy[x].direction[i].blue;
1480 if (image->colorspace == CMYKColorspace)
1481 sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1482 density_xy[x].direction[i].black;
1483 if (image->alpha_trait != UndefinedPixelTrait)
1484 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1485 density_xy[x].direction[i].alpha;
1486 /*
1487 Difference entropy.
1488 */
1489 channel_features[RedPixelChannel].difference_entropy[i]-=
1490 density_xy[x].direction[i].red*
1491 MagickLog10(density_xy[x].direction[i].red);
1492 channel_features[GreenPixelChannel].difference_entropy[i]-=
1493 density_xy[x].direction[i].green*
1494 MagickLog10(density_xy[x].direction[i].green);
1495 channel_features[BluePixelChannel].difference_entropy[i]-=
1496 density_xy[x].direction[i].blue*
1497 MagickLog10(density_xy[x].direction[i].blue);
1498 if (image->colorspace == CMYKColorspace)
1499 channel_features[BlackPixelChannel].difference_entropy[i]-=
1500 density_xy[x].direction[i].black*
1501 MagickLog10(density_xy[x].direction[i].black);
1502 if (image->alpha_trait != UndefinedPixelTrait)
1503 channel_features[AlphaPixelChannel].difference_entropy[i]-=
1504 density_xy[x].direction[i].alpha*
1505 MagickLog10(density_xy[x].direction[i].alpha);
1506 /*
1507 Information Measures of Correlation.
1508 */
1509 entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1510 MagickLog10(density_x[x].direction[i].red));
1511 entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1512 MagickLog10(density_x[x].direction[i].green));
1513 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1514 MagickLog10(density_x[x].direction[i].blue));
1515 if (image->colorspace == CMYKColorspace)
1516 entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1517 MagickLog10(density_x[x].direction[i].black));
1518 if (image->alpha_trait != UndefinedPixelTrait)
1519 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1520 MagickLog10(density_x[x].direction[i].alpha));
1521 entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1522 MagickLog10(density_y[x].direction[i].red));
1523 entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1524 MagickLog10(density_y[x].direction[i].green));
1525 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1526 MagickLog10(density_y[x].direction[i].blue));
1527 if (image->colorspace == CMYKColorspace)
1528 entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1529 MagickLog10(density_y[x].direction[i].black));
1530 if (image->alpha_trait != UndefinedPixelTrait)
1531 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1532 MagickLog10(density_y[x].direction[i].alpha));
1533 }
1534 /*
1535 Difference variance.
1536 */
1537 channel_features[RedPixelChannel].difference_variance[i]=
1538 (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1539 (variance.direction[i].red*variance.direction[i].red))/
1540 ((double) number_grays*number_grays*number_grays*number_grays);
1541 channel_features[GreenPixelChannel].difference_variance[i]=
1542 (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1543 (variance.direction[i].green*variance.direction[i].green))/
1544 ((double) number_grays*number_grays*number_grays*number_grays);
1545 channel_features[BluePixelChannel].difference_variance[i]=
1546 (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1547 (variance.direction[i].blue*variance.direction[i].blue))/
1548 ((double) number_grays*number_grays*number_grays*number_grays);
1549 if (image->colorspace == CMYKColorspace)
1550 channel_features[BlackPixelChannel].difference_variance[i]=
1551 (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1552 (variance.direction[i].black*variance.direction[i].black))/
1553 ((double) number_grays*number_grays*number_grays*number_grays);
1554 if (image->alpha_trait != UndefinedPixelTrait)
1555 channel_features[AlphaPixelChannel].difference_variance[i]=
1556 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1557 (variance.direction[i].alpha*variance.direction[i].alpha))/
1558 ((double) number_grays*number_grays*number_grays*number_grays);
1559 /*
1560 Information Measures of Correlation.
1561 */
1562 channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1563 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1564 (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1565 entropy_x.direction[i].red : entropy_y.direction[i].red);
1566 channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1567 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1568 (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1569 entropy_x.direction[i].green : entropy_y.direction[i].green);
1570 channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1571 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1572 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1573 entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1574 if (image->colorspace == CMYKColorspace)
1575 channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1576 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1577 (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1578 entropy_x.direction[i].black : entropy_y.direction[i].black);
1579 if (image->alpha_trait != UndefinedPixelTrait)
1580 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1581 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1582 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1583 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1584 channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1585 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].red-
1586 entropy_xy.direction[i].red)))));
1587 channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1588 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].green-
1589 entropy_xy.direction[i].green)))));
1590 channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1591 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].blue-
1592 entropy_xy.direction[i].blue)))));
1593 if (image->colorspace == CMYKColorspace)
1594 channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1595 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].black-
1596 entropy_xy.direction[i].black)))));
1597 if (image->alpha_trait != UndefinedPixelTrait)
1598 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1599 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha-
1600 entropy_xy.direction[i].alpha)))));
1601 }
1602 /*
1603 Compute more texture features.
1604 */
1605 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1606 #pragma omp parallel for schedule(static) shared(status) \
1607 magick_number_threads(image,image,number_grays,1)
1608 #endif
1609 for (i=0; i < 4; i++)
1610 {
1611 ssize_t
1612 z;
1613
1614 for (z=0; z < (ssize_t) number_grays; z++)
1615 {
1616 ssize_t
1617 y;
1618
1619 ChannelStatistics
1620 pixel;
1621
1622 (void) memset(&pixel,0,sizeof(pixel));
1623 for (y=0; y < (ssize_t) number_grays; y++)
1624 {
1625 ssize_t
1626 x;
1627
1628 for (x=0; x < (ssize_t) number_grays; x++)
1629 {
1630 /*
1631 Contrast: amount of local variations present in an image.
1632 */
1633 if (((y-x) == z) || ((x-y) == z))
1634 {
1635 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1636 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1637 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1638 if (image->colorspace == CMYKColorspace)
1639 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1640 if (image->alpha_trait != UndefinedPixelTrait)
1641 pixel.direction[i].alpha+=
1642 cooccurrence[x][y].direction[i].alpha;
1643 }
1644 /*
1645 Maximum Correlation Coefficient.
1646 */
1647 if ((fabs(density_x[z].direction[i].red) > MagickEpsilon) &&
1648 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1649 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1650 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1651 density_y[x].direction[i].red;
1652 if ((fabs(density_x[z].direction[i].green) > MagickEpsilon) &&
1653 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1654 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1655 cooccurrence[y][x].direction[i].green/
1656 density_x[z].direction[i].green/density_y[x].direction[i].red;
1657 if ((fabs(density_x[z].direction[i].blue) > MagickEpsilon) &&
1658 (fabs(density_y[x].direction[i].blue) > MagickEpsilon))
1659 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1660 cooccurrence[y][x].direction[i].blue/
1661 density_x[z].direction[i].blue/density_y[x].direction[i].blue;
1662 if (image->colorspace == CMYKColorspace)
1663 if ((fabs(density_x[z].direction[i].black) > MagickEpsilon) &&
1664 (fabs(density_y[x].direction[i].black) > MagickEpsilon))
1665 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1666 cooccurrence[y][x].direction[i].black/
1667 density_x[z].direction[i].black/density_y[x].direction[i].black;
1668 if (image->alpha_trait != UndefinedPixelTrait)
1669 if ((fabs(density_x[z].direction[i].alpha) > MagickEpsilon) &&
1670 (fabs(density_y[x].direction[i].alpha) > MagickEpsilon))
1671 Q[z][y].direction[i].alpha+=
1672 cooccurrence[z][x].direction[i].alpha*
1673 cooccurrence[y][x].direction[i].alpha/
1674 density_x[z].direction[i].alpha/
1675 density_y[x].direction[i].alpha;
1676 }
1677 }
1678 channel_features[RedPixelChannel].contrast[i]+=z*z*
1679 pixel.direction[i].red;
1680 channel_features[GreenPixelChannel].contrast[i]+=z*z*
1681 pixel.direction[i].green;
1682 channel_features[BluePixelChannel].contrast[i]+=z*z*
1683 pixel.direction[i].blue;
1684 if (image->colorspace == CMYKColorspace)
1685 channel_features[BlackPixelChannel].contrast[i]+=z*z*
1686 pixel.direction[i].black;
1687 if (image->alpha_trait != UndefinedPixelTrait)
1688 channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1689 pixel.direction[i].alpha;
1690 }
1691 /*
1692 Maximum Correlation Coefficient.
1693 Future: return second largest eigenvalue of Q.
1694 */
1695 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1696 sqrt((double) -1.0);
1697 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1698 sqrt((double) -1.0);
1699 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1700 sqrt((double) -1.0);
1701 if (image->colorspace == CMYKColorspace)
1702 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1703 sqrt((double) -1.0);
1704 if (image->alpha_trait != UndefinedPixelTrait)
1705 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1706 sqrt((double) -1.0);
1707 }
1708 /*
1709 Relinquish resources.
1710 */
1711 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1712 for (i=0; i < (ssize_t) number_grays; i++)
1713 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1714 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1715 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1716 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1717 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1718 for (i=0; i < (ssize_t) number_grays; i++)
1719 cooccurrence[i]=(ChannelStatistics *)
1720 RelinquishMagickMemory(cooccurrence[i]);
1721 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1722 return(channel_features);
1723 }
1724
1725 /*
1726 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1727 % %
1728 % %
1729 % %
1730 % H o u g h L i n e I m a g e %
1731 % %
1732 % %
1733 % %
1734 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1735 %
1736 % Use HoughLineImage() in conjunction with any binary edge extracted image (we
1737 % recommand Canny) to identify lines in the image. The algorithm accumulates
1738 % counts for every white pixel for every possible orientation (for angles from
1739 % 0 to 179 in 1 degree increments) and distance from the center of the image to
1740 % the corner (in 1 px increments) and stores the counts in an accumulator
1741 % matrix of angle vs distance. The size of the accumulator is 180x(diagonal/2).
1742 % Next it searches this space for peaks in counts and converts the locations
1743 % of the peaks to slope and intercept in the normal x,y input image space. Use
1744 % the slope/intercepts to find the endpoints clipped to the bounds of the
1745 % image. The lines are then drawn. The counts are a measure of the length of
1746 % the lines.
1747 %
1748 % The format of the HoughLineImage method is:
1749 %
1750 % Image *HoughLineImage(const Image *image,const size_t width,
1751 % const size_t height,const size_t threshold,ExceptionInfo *exception)
1752 %
1753 % A description of each parameter follows:
1754 %
1755 % o image: the image.
1756 %
1757 % o width, height: find line pairs as local maxima in this neighborhood.
1758 %
1759 % o threshold: the line count threshold.
1760 %
1761 % o exception: return any errors or warnings in this structure.
1762 %
1763 */
1764
MagickRound(double x)1765 static inline double MagickRound(double x)
1766 {
1767 /*
1768 Round the fraction to nearest integer.
1769 */
1770 if ((x-floor(x)) < (ceil(x)-x))
1771 return(floor(x));
1772 return(ceil(x));
1773 }
1774
RenderHoughLines(const ImageInfo * image_info,const size_t columns,const size_t rows,ExceptionInfo * exception)1775 static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns,
1776 const size_t rows,ExceptionInfo *exception)
1777 {
1778 #define BoundingBox "viewbox"
1779
1780 DrawInfo
1781 *draw_info;
1782
1783 Image
1784 *image;
1785
1786 MagickBooleanType
1787 status;
1788
1789 /*
1790 Open image.
1791 */
1792 image=AcquireImage(image_info,exception);
1793 status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
1794 if (status == MagickFalse)
1795 {
1796 image=DestroyImageList(image);
1797 return((Image *) NULL);
1798 }
1799 image->columns=columns;
1800 image->rows=rows;
1801 draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL);
1802 draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/
1803 DefaultResolution;
1804 draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/
1805 DefaultResolution;
1806 image->columns=(size_t) (draw_info->affine.sx*image->columns);
1807 image->rows=(size_t) (draw_info->affine.sy*image->rows);
1808 status=SetImageExtent(image,image->columns,image->rows,exception);
1809 if (status == MagickFalse)
1810 return(DestroyImageList(image));
1811 if (SetImageBackgroundColor(image,exception) == MagickFalse)
1812 {
1813 image=DestroyImageList(image);
1814 return((Image *) NULL);
1815 }
1816 /*
1817 Render drawing.
1818 */
1819 if (GetBlobStreamData(image) == (unsigned char *) NULL)
1820 draw_info->primitive=FileToString(image->filename,~0UL,exception);
1821 else
1822 {
1823 draw_info->primitive=(char *) AcquireQuantumMemory(1,(size_t)
1824 GetBlobSize(image)+1);
1825 if (draw_info->primitive != (char *) NULL)
1826 {
1827 (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
1828 (size_t) GetBlobSize(image));
1829 draw_info->primitive[GetBlobSize(image)]='\0';
1830 }
1831 }
1832 (void) DrawImage(image,draw_info,exception);
1833 draw_info=DestroyDrawInfo(draw_info);
1834 (void) CloseBlob(image);
1835 return(GetFirstImageInList(image));
1836 }
1837
HoughLineImage(const Image * image,const size_t width,const size_t height,const size_t threshold,ExceptionInfo * exception)1838 MagickExport Image *HoughLineImage(const Image *image,const size_t width,
1839 const size_t height,const size_t threshold,ExceptionInfo *exception)
1840 {
1841 #define HoughLineImageTag "HoughLine/Image"
1842
1843 CacheView
1844 *image_view;
1845
1846 char
1847 message[MagickPathExtent],
1848 path[MagickPathExtent];
1849
1850 const char
1851 *artifact;
1852
1853 double
1854 hough_height;
1855
1856 Image
1857 *lines_image = NULL;
1858
1859 ImageInfo
1860 *image_info;
1861
1862 int
1863 file;
1864
1865 MagickBooleanType
1866 status;
1867
1868 MagickOffsetType
1869 progress;
1870
1871 MatrixInfo
1872 *accumulator;
1873
1874 PointInfo
1875 center;
1876
1877 ssize_t
1878 y;
1879
1880 size_t
1881 accumulator_height,
1882 accumulator_width,
1883 line_count;
1884
1885 /*
1886 Create the accumulator.
1887 */
1888 assert(image != (const Image *) NULL);
1889 assert(image->signature == MagickCoreSignature);
1890 if (image->debug != MagickFalse)
1891 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1892 assert(exception != (ExceptionInfo *) NULL);
1893 assert(exception->signature == MagickCoreSignature);
1894 accumulator_width=180;
1895 hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1896 image->rows : image->columns))/2.0);
1897 accumulator_height=(size_t) (2.0*hough_height);
1898 accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1899 sizeof(double),exception);
1900 if (accumulator == (MatrixInfo *) NULL)
1901 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1902 if (NullMatrix(accumulator) == MagickFalse)
1903 {
1904 accumulator=DestroyMatrixInfo(accumulator);
1905 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1906 }
1907 /*
1908 Populate the accumulator.
1909 */
1910 status=MagickTrue;
1911 progress=0;
1912 center.x=(double) image->columns/2.0;
1913 center.y=(double) image->rows/2.0;
1914 image_view=AcquireVirtualCacheView(image,exception);
1915 for (y=0; y < (ssize_t) image->rows; y++)
1916 {
1917 const Quantum
1918 *magick_restrict p;
1919
1920 ssize_t
1921 x;
1922
1923 if (status == MagickFalse)
1924 continue;
1925 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1926 if (p == (Quantum *) NULL)
1927 {
1928 status=MagickFalse;
1929 continue;
1930 }
1931 for (x=0; x < (ssize_t) image->columns; x++)
1932 {
1933 if (GetPixelIntensity(image,p) > (QuantumRange/2.0))
1934 {
1935 ssize_t
1936 i;
1937
1938 for (i=0; i < 180; i++)
1939 {
1940 double
1941 count,
1942 radius;
1943
1944 radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
1945 (((double) y-center.y)*sin(DegreesToRadians((double) i)));
1946 (void) GetMatrixElement(accumulator,i,(ssize_t)
1947 MagickRound(radius+hough_height),&count);
1948 count++;
1949 (void) SetMatrixElement(accumulator,i,(ssize_t)
1950 MagickRound(radius+hough_height),&count);
1951 }
1952 }
1953 p+=GetPixelChannels(image);
1954 }
1955 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1956 {
1957 MagickBooleanType
1958 proceed;
1959
1960 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1961 #pragma omp atomic
1962 #endif
1963 progress++;
1964 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
1965 if (proceed == MagickFalse)
1966 status=MagickFalse;
1967 }
1968 }
1969 image_view=DestroyCacheView(image_view);
1970 if (status == MagickFalse)
1971 {
1972 accumulator=DestroyMatrixInfo(accumulator);
1973 return((Image *) NULL);
1974 }
1975 /*
1976 Generate line segments from accumulator.
1977 */
1978 file=AcquireUniqueFileResource(path);
1979 if (file == -1)
1980 {
1981 accumulator=DestroyMatrixInfo(accumulator);
1982 return((Image *) NULL);
1983 }
1984 (void) FormatLocaleString(message,MagickPathExtent,
1985 "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
1986 (double) height,(double) threshold);
1987 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1988 status=MagickFalse;
1989 (void) FormatLocaleString(message,MagickPathExtent,
1990 "viewbox 0 0 %.20g %.20g\n",(double) image->columns,(double) image->rows);
1991 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1992 status=MagickFalse;
1993 (void) FormatLocaleString(message,MagickPathExtent,
1994 "# x1,y1 x2,y2 # count angle distance\n");
1995 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1996 status=MagickFalse;
1997 line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
1998 if (threshold != 0)
1999 line_count=threshold;
2000 for (y=0; y < (ssize_t) accumulator_height; y++)
2001 {
2002 ssize_t
2003 x;
2004
2005 for (x=0; x < (ssize_t) accumulator_width; x++)
2006 {
2007 double
2008 count;
2009
2010 (void) GetMatrixElement(accumulator,x,y,&count);
2011 if (count >= (double) line_count)
2012 {
2013 double
2014 maxima;
2015
2016 SegmentInfo
2017 line;
2018
2019 ssize_t
2020 v;
2021
2022 /*
2023 Is point a local maxima?
2024 */
2025 maxima=count;
2026 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2027 {
2028 ssize_t
2029 u;
2030
2031 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2032 {
2033 if ((u != 0) || (v !=0))
2034 {
2035 (void) GetMatrixElement(accumulator,x+u,y+v,&count);
2036 if (count > maxima)
2037 {
2038 maxima=count;
2039 break;
2040 }
2041 }
2042 }
2043 if (u < (ssize_t) (width/2))
2044 break;
2045 }
2046 (void) GetMatrixElement(accumulator,x,y,&count);
2047 if (maxima > count)
2048 continue;
2049 if ((x >= 45) && (x <= 135))
2050 {
2051 /*
2052 y = (r-x cos(t))/sin(t)
2053 */
2054 line.x1=0.0;
2055 line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
2056 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2057 sin(DegreesToRadians((double) x))+(image->rows/2.0);
2058 line.x2=(double) image->columns;
2059 line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
2060 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2061 sin(DegreesToRadians((double) x))+(image->rows/2.0);
2062 }
2063 else
2064 {
2065 /*
2066 x = (r-y cos(t))/sin(t)
2067 */
2068 line.y1=0.0;
2069 line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2070 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2071 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2072 line.y2=(double) image->rows;
2073 line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2074 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2075 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2076 }
2077 (void) FormatLocaleString(message,MagickPathExtent,
2078 "line %g,%g %g,%g # %g %g %g\n",line.x1,line.y1,line.x2,line.y2,
2079 maxima,(double) x,(double) y);
2080 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2081 status=MagickFalse;
2082 }
2083 }
2084 }
2085 (void) close(file);
2086 /*
2087 Render lines to image canvas.
2088 */
2089 image_info=AcquireImageInfo();
2090 image_info->background_color=image->background_color;
2091 (void) FormatLocaleString(image_info->filename,MagickPathExtent,"%s",path);
2092 artifact=GetImageArtifact(image,"background");
2093 if (artifact != (const char *) NULL)
2094 (void) SetImageOption(image_info,"background",artifact);
2095 artifact=GetImageArtifact(image,"fill");
2096 if (artifact != (const char *) NULL)
2097 (void) SetImageOption(image_info,"fill",artifact);
2098 artifact=GetImageArtifact(image,"stroke");
2099 if (artifact != (const char *) NULL)
2100 (void) SetImageOption(image_info,"stroke",artifact);
2101 artifact=GetImageArtifact(image,"strokewidth");
2102 if (artifact != (const char *) NULL)
2103 (void) SetImageOption(image_info,"strokewidth",artifact);
2104 lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception);
2105 artifact=GetImageArtifact(image,"hough-lines:accumulator");
2106 if ((lines_image != (Image *) NULL) &&
2107 (IsStringTrue(artifact) != MagickFalse))
2108 {
2109 Image
2110 *accumulator_image;
2111
2112 accumulator_image=MatrixToImage(accumulator,exception);
2113 if (accumulator_image != (Image *) NULL)
2114 AppendImageToList(&lines_image,accumulator_image);
2115 }
2116 /*
2117 Free resources.
2118 */
2119 accumulator=DestroyMatrixInfo(accumulator);
2120 image_info=DestroyImageInfo(image_info);
2121 (void) RelinquishUniqueFileResource(path);
2122 return(GetFirstImageInList(lines_image));
2123 }
2124
2125 /*
2126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2127 % %
2128 % %
2129 % %
2130 % M e a n S h i f t I m a g e %
2131 % %
2132 % %
2133 % %
2134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2135 %
2136 % MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
2137 % each pixel, it visits all the pixels in the neighborhood specified by
2138 % the window centered at the pixel and excludes those that are outside the
2139 % radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
2140 % that are within the specified color distance from the current mean, and
2141 % computes a new x,y centroid from those coordinates and a new mean. This new
2142 % x,y centroid is used as the center for a new window. This process iterates
2143 % until it converges and the final mean is replaces the (original window
2144 % center) pixel value. It repeats this process for the next pixel, etc.,
2145 % until it processes all pixels in the image. Results are typically better with
2146 % colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
2147 %
2148 % The format of the MeanShiftImage method is:
2149 %
2150 % Image *MeanShiftImage(const Image *image,const size_t width,
2151 % const size_t height,const double color_distance,
2152 % ExceptionInfo *exception)
2153 %
2154 % A description of each parameter follows:
2155 %
2156 % o image: the image.
2157 %
2158 % o width, height: find pixels in this neighborhood.
2159 %
2160 % o color_distance: the color distance.
2161 %
2162 % o exception: return any errors or warnings in this structure.
2163 %
2164 */
MeanShiftImage(const Image * image,const size_t width,const size_t height,const double color_distance,ExceptionInfo * exception)2165 MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
2166 const size_t height,const double color_distance,ExceptionInfo *exception)
2167 {
2168 #define MaxMeanShiftIterations 100
2169 #define MeanShiftImageTag "MeanShift/Image"
2170
2171 CacheView
2172 *image_view,
2173 *mean_view,
2174 *pixel_view;
2175
2176 Image
2177 *mean_image;
2178
2179 MagickBooleanType
2180 status;
2181
2182 MagickOffsetType
2183 progress;
2184
2185 ssize_t
2186 y;
2187
2188 assert(image != (const Image *) NULL);
2189 assert(image->signature == MagickCoreSignature);
2190 if (image->debug != MagickFalse)
2191 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2192 assert(exception != (ExceptionInfo *) NULL);
2193 assert(exception->signature == MagickCoreSignature);
2194 mean_image=CloneImage(image,0,0,MagickTrue,exception);
2195 if (mean_image == (Image *) NULL)
2196 return((Image *) NULL);
2197 if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
2198 {
2199 mean_image=DestroyImage(mean_image);
2200 return((Image *) NULL);
2201 }
2202 status=MagickTrue;
2203 progress=0;
2204 image_view=AcquireVirtualCacheView(image,exception);
2205 pixel_view=AcquireVirtualCacheView(image,exception);
2206 mean_view=AcquireAuthenticCacheView(mean_image,exception);
2207 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2208 #pragma omp parallel for schedule(static) shared(status,progress) \
2209 magick_number_threads(mean_image,mean_image,mean_image->rows,1)
2210 #endif
2211 for (y=0; y < (ssize_t) mean_image->rows; y++)
2212 {
2213 const Quantum
2214 *magick_restrict p;
2215
2216 Quantum
2217 *magick_restrict q;
2218
2219 ssize_t
2220 x;
2221
2222 if (status == MagickFalse)
2223 continue;
2224 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2225 q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2226 exception);
2227 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2228 {
2229 status=MagickFalse;
2230 continue;
2231 }
2232 for (x=0; x < (ssize_t) mean_image->columns; x++)
2233 {
2234 PixelInfo
2235 mean_pixel,
2236 previous_pixel;
2237
2238 PointInfo
2239 mean_location,
2240 previous_location;
2241
2242 ssize_t
2243 i;
2244
2245 GetPixelInfo(image,&mean_pixel);
2246 GetPixelInfoPixel(image,p,&mean_pixel);
2247 mean_location.x=(double) x;
2248 mean_location.y=(double) y;
2249 for (i=0; i < MaxMeanShiftIterations; i++)
2250 {
2251 double
2252 distance,
2253 gamma;
2254
2255 PixelInfo
2256 sum_pixel;
2257
2258 PointInfo
2259 sum_location;
2260
2261 ssize_t
2262 count,
2263 v;
2264
2265 sum_location.x=0.0;
2266 sum_location.y=0.0;
2267 GetPixelInfo(image,&sum_pixel);
2268 previous_location=mean_location;
2269 previous_pixel=mean_pixel;
2270 count=0;
2271 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2272 {
2273 ssize_t
2274 u;
2275
2276 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2277 {
2278 if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2279 {
2280 PixelInfo
2281 pixel;
2282
2283 status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
2284 MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2285 mean_location.y+v),&pixel,exception);
2286 distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
2287 (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
2288 (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
2289 if (distance <= (color_distance*color_distance))
2290 {
2291 sum_location.x+=mean_location.x+u;
2292 sum_location.y+=mean_location.y+v;
2293 sum_pixel.red+=pixel.red;
2294 sum_pixel.green+=pixel.green;
2295 sum_pixel.blue+=pixel.blue;
2296 sum_pixel.alpha+=pixel.alpha;
2297 count++;
2298 }
2299 }
2300 }
2301 }
2302 gamma=PerceptibleReciprocal(count);
2303 mean_location.x=gamma*sum_location.x;
2304 mean_location.y=gamma*sum_location.y;
2305 mean_pixel.red=gamma*sum_pixel.red;
2306 mean_pixel.green=gamma*sum_pixel.green;
2307 mean_pixel.blue=gamma*sum_pixel.blue;
2308 mean_pixel.alpha=gamma*sum_pixel.alpha;
2309 distance=(mean_location.x-previous_location.x)*
2310 (mean_location.x-previous_location.x)+
2311 (mean_location.y-previous_location.y)*
2312 (mean_location.y-previous_location.y)+
2313 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2314 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2315 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2316 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2317 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2318 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2319 if (distance <= 3.0)
2320 break;
2321 }
2322 SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
2323 SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
2324 SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
2325 SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
2326 p+=GetPixelChannels(image);
2327 q+=GetPixelChannels(mean_image);
2328 }
2329 if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2330 status=MagickFalse;
2331 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2332 {
2333 MagickBooleanType
2334 proceed;
2335
2336 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2337 #pragma omp atomic
2338 #endif
2339 progress++;
2340 proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
2341 if (proceed == MagickFalse)
2342 status=MagickFalse;
2343 }
2344 }
2345 mean_view=DestroyCacheView(mean_view);
2346 pixel_view=DestroyCacheView(pixel_view);
2347 image_view=DestroyCacheView(image_view);
2348 return(mean_image);
2349 }
2350