1/* Various operators as functions. 2 */ 3 4logical_and a b = a && b; 5logical_or a b = a || b; 6bitwise_and a b = a & b; 7bitwise_or a b = a | b; 8eor a b = a ^ b; 9left_shift a b = a << b; 10right_shift a b = a >> b; 11not a = !a; 12 13less a b = a < b; 14more a b = a > b; 15less_equal a b = a <= b; 16more_equal a b = a >= b; 17equal a b = a == b; 18not_equal a b = a != b; 19pointer_equal a b = a === b; 20not_pointer_equal a b = a !== b; 21 22add a b = a + b; 23subtract a b = a - b; 24multiply a b = a * b; 25divide a b = a / b; 26idivide a b = (int) ((int) a / (int) b); 27power a b = a ** b; 28square x = x * x; 29remainder a b = a % b; 30 31cons a b = a : b; 32dot a b = a . ( b ); 33join a b = a ++ b; 34subscript a b = a ? b; 35 36generate s n f = [s, n .. f]; 37comma r i = (r, i); 38 39compose f g = f @ g; 40 41// our only trinary operator is actually a binary operator 42if_then_else a x = if a then x?0 else x?1; 43 44cast_unsigned_char x = (unsigned char) x; 45cast_signed_char x = (signed char) x; 46cast_unsigned_short x = (unsigned short) x; 47cast_signed_short x = (signed short) x; 48cast_unsigned_int x = (unsigned int) x; 49cast_signed_int x = (signed int) x; 50cast_float x = (float) x; 51cast_double x = (double) x; 52cast_complex x = (complex) x; 53cast_double_complex x = (double complex) x; 54 55unary_minus x = -x; 56negate x = !x; 57complement x = ~x; 58unary_plus x = +x; 59 60// the function we call for "a -> v" expressions 61mksvpair s v 62 = [s, v], is_string s 63 = error "not str on lhs of ->"; 64 65// the vector ops ... im is an image, vec is a real_list 66vec op_name im vec 67 = im_lintra_vec ones im vec, 68 op_name == "add" || op_name == "add'" 69 = im_lintra_vec ones (-1 * im) vec, 70 op_name == "subtract'" 71 = im_lintra_vec ones im inv, 72 op_name == "subtract" 73 = im_lintra_vec vec im zeros, 74 op_name == "multiply" || op_name == "multiply'" 75 = im_lintra_vec vec (1 / im) zeros, 76 op_name == "divide'" 77 = im_lintra_vec recip im zeros, 78 op_name == "divide" 79 = im_expntra_vec im vec, 80 op_name == "power'" 81 = im_powtra_vec im vec, 82 op_name == "power" 83 = im_remainderconst_vec im vec, 84 op_name == "remainder" 85 = im_andimage_vec im vec, 86 op_name == "bitwise_and" || op_name == "bitwise_and'" 87 = im_orimage_vec im vec, 88 op_name == "bitwise_or" || op_name == "bitwise_or'" 89 = im_eorimage_vec im vec, 90 op_name == "eor" || op_name == "eor'" 91 = im_equal_vec im vec, 92 op_name == "equal" || op_name == "equal'" 93 = im_notequal_vec im vec, 94 op_name == "not_equal" || op_name == "not_equal'" 95 = im_less_vec im vec, 96 op_name == "less" 97 = im_moreeq_vec im vec, 98 op_name == "less'" 99 = im_lesseq_vec im vec, 100 op_name == "less_equal" 101 = im_more_vec im vec, 102 op_name == "less_equal'" 103 = error "unimplemented vector operation" 104{ 105 zeros = replicate (len vec) 0; 106 ones = replicate (len vec) 1; 107 recip = map (divide 1) vec; 108 inv = map (multiply (-1)) vec; 109} 110 111// make a name value pair 112mknvpair n v 113 = [n, v], is_string n 114 = error "not [char] on LHS of =>"; 115 116/* Macbeth chart patch names. 117 */ 118macbeth_names = [ 119 "Dark skin", 120 "Light skin", 121 "Blue sky", 122 "Foliage", 123 "Blue flower", 124 "Bluish green", 125 "Orange", 126 "Purplish blue", 127 "Moderate red", 128 "Purple", 129 "Yellow green", 130 "Orange yellow", 131 "Blue", 132 "Green", 133 "Red", 134 "Yellow", 135 "Magenta", 136 "Cyan", 137 "White (density 0.05)", 138 "Neutral 8 (density 0.23)", 139 "Neutral 6.5 (density 0.44)", 140 "Neutral 5 (density 0.70)", 141 "Neutral 3.5 (density 1.05)", 142 "Black (density 1.50)" 143]; 144 145bandsplit x 146 = oo_unary_function bandsplit_op x, is_class x 147 = map (subscript x) [0 .. bands - 1], is_image x 148 = error (_ "bad arguments to " ++ "bandsplit") 149{ 150 bands = get_header "Bands" x; 151 bandsplit_op = Operator "bandsplit" (map Image @ bandsplit) 152 Operator_type.COMPOUND false; 153} 154 155bandjoin l 156 = wrapper joined, has_wrapper 157 = joined, is_listof has_image l 158 = error (_ "bad arguments to " ++ "bandjoin") 159{ 160 has_wrapper = has_member_list (has_member "Image") l; 161 wrapper = get_member_list (has_member "Image") (get_member "Image") l; 162 joined = im_gbandjoin (map get_image l); 163} 164 165mean x 166 = oo_unary_function mean_op x, is_class x 167 = im_avg x, is_image x 168 = mean_list x, is_list x 169 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "mean") 170{ 171 mean_op = Operator "mean" mean Operator_type.COMPOUND false; 172 173 mean_list l = sum l / size l; 174 175 // number of elements in some sort of nested-list thing 176 size l 177 = foldr acc 0 l 178 { 179 acc x total 180 = total + size x, is_list x 181 = total + 1; 182 } 183 184 // add elements in a nested-list thing 185 sum l 186 = foldr acc 0 l 187 { 188 acc x total 189 = total + sum x, is_list x 190 = total + x; 191 } 192} 193 194meang x 195 = (appl (power e) @ mean @ appl log) x 196{ 197 appl fn x 198 = map fn x, is_list x 199 = fn x; 200} 201 202// zero-excluding mean 203meanze x 204 = oo_unary_function meanze_op x, is_class x 205 = meanze_image x, is_image x 206 = meanze_list x, is_list x 207 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "meanze") 208{ 209 meanze_op = Operator "meanze" meanze Operator_type.COMPOUND false; 210 211 meanze_list l = sum l / size l; 212 213 // number of non-zero elements in some sort of nested-list thing 214 size l 215 = foldr acc 0 l 216 { 217 acc x total 218 = total + size x, is_list x 219 = total + 1, x != 0; 220 = total; 221 } 222 223 // add elements in a nested-list thing 224 sum l 225 = foldr acc 0 l 226 { 227 acc x total 228 = total + sum x, is_list x 229 = total + x; 230 } 231 232 // image mean 233 meanze_image i 234 = sum / N 235 { 236 w = get_width i; 237 h = get_height i; 238 b = get_bands i; 239 240 st = stats i; 241 sum = st.value?0?2; 242 243 // find non-zero pixels (not zero in all bands) 244 zp = im_notequal_vec i (replicate b 0); 245 246 // number of non-zero pixels 247 N = b * (mean zp * w * h) / 255; 248 } 249} 250 251deviation x 252 = oo_unary_function deviation_op x, is_class x 253 = im_deviate x, is_image x 254 = deviation_list x, is_real_list x || is_matrix x 255 = error (_ "bad arguments to " ++ "deviation") 256{ 257 deviation_op = Operator "deviation" 258 deviation Operator_type.COMPOUND false; 259 260 deviation_list l 261 = (abs (s2 - (s * s / n)) / (n - 1)) ** 0.5 262 { 263 [n, s, s2] = sum_sum2_list l; 264 } 265 266 // return n, sum, sum of squares for a list of reals 267 sum_sum2_list x 268 = foldr accumulate [0, 0, 0] x 269 { 270 accumulate x sofar 271 = [n + 1, x + s, x * x + s2], is_real x 272 = [n + n', s + s', s2 + s2'], is_list x 273 = error "sum_sum2_list: not real or [real]" 274 { 275 [n, s, s2] = sofar; 276 [n', s', s2'] = sum_sum2_list x; 277 } 278 } 279} 280 281deviationze x 282 = oo_unary_function deviationze_op x, is_class x 283 = deviationze_image x, is_image x 284 = deviationze_list x, is_real_list x || is_matrix x 285 = error (_ "bad arguments to " ++ "deviationze") 286{ 287 deviationze_op = Operator "deviationze" 288 deviationze Operator_type.COMPOUND false; 289 290 deviationze_list l 291 = (abs (s2 - (s * s / n)) / (n - 1)) ** 0.5 292 { 293 [n, s, s2] = sum_sum2_list l; 294 } 295 296 // return number of non-zero elements, sum, sum of squares for a list of 297 // reals 298 sum_sum2_list x 299 = foldr accumulate [0, 0, 0] x 300 { 301 accumulate x sofar 302 = sofar, is_real x && x == 0 303 = [n + 1, x + s, x * x + s2], is_real x 304 = [n + n', s + s', s2 + s2'], is_list x 305 = error "sum_sum2_list: not real or [real]" 306 { 307 [n, s, s2] = sofar; 308 [n', s', s2'] = sum_sum2_list x; 309 } 310 } 311 312 deviationze_image i 313 = ((sum2 - sum * sum / N) / (N - 1)) ** 0.5 314 { 315 w = get_width i; 316 h = get_height i; 317 b = get_bands i; 318 319 st = stats i; 320 sum = st.value?0?2; 321 sum2 = st.value?0?3; 322 323 // find non-zero pixels (not zero in all bands) 324 zp = im_notequal_vec i (replicate b 0); 325 326 // number of non-zero pixels 327 N = b * (mean zp * w * h) / 255; 328 } 329} 330 331// find the centre of gravity of a histogram 332gravity x 333 = oo_unary_function gravity_op x, is_class x 334 = im_hist_gravity x, is_hist x 335 = gravity_list x, is_list x 336 = error (_ "bad arguments to " ++ "gravity") 337{ 338 gravity_op = Operator "gravity" gravity Operator_type.COMPOUND false; 339 340 // centre of gravity of a histogram... use the histogram to weight an 341 // identity, then sum, then find the mean element 342 im_hist_gravity h 343 = m 344 { 345 // make horizontal 346 h' 347 = rot270 h, get_width h == 1 348 = h, get_height h == 1 349 = error "width or height not 1"; 350 351 // number of elements 352 w = get_width h'; 353 354 // matching identity 355 i 356 = im_identity_ushort 1 w, w <= 2 ** 16 - 1 357 = make_xy w 1 ? 0; 358 359 // weight identity and sum 360 s = mean (i * h') * w; 361 362 // sum of original histogram 363 s' = mean h * w; 364 365 // weighted mean 366 m = s / s'; 367 } 368 369 gravity_list l 370 = m 371 { 372 w = len l; 373 374 // matching identity 375 i = [0, 1 .. w - 1]; 376 377 // weight identity and sum 378 s = sum (map2 multiply i l); 379 380 // sum of original histogram 381 s' = sum l; 382 383 // weighted mean 384 m = s / s'; 385 } 386} 387 388project x 389 = oo_unary_function project_op x, is_class x 390 = im_project x, is_image x 391 = error (_ "bad arguments to " ++ "project") 392{ 393 project_op = Operator "project" project Operator_type.COMPOUND false; 394} 395 396abs x 397 = oo_unary_function abs_op x, is_class x 398 = im_abs x, is_image x 399 = abs_cmplx x, is_complex x 400 = abs_num x, is_real x 401 = abs_list x, is_real_list x 402 = abs_list (map abs_list x), is_matrix x 403 = error (_ "bad arguments to " ++ "abs") 404{ 405 abs_op = Operator "abs" abs Operator_type.COMPOUND false; 406 407 abs_list l = (sum (map square l)) ** 0.5; 408 409 abs_num n 410 = n, n >= 0 411 = -n; 412 413 abs_cmplx c = ((re c)**2 + (im c)**2) ** 0.5; 414} 415 416copy x 417 = oo_unary_function copy_op x, is_class x 418 = im_copy x, is_image x 419 = x 420{ 421 copy_op = Operator "copy" copy Operator_type.COMPOUND_REWRAP false; 422} 423 424// like abs, but treat pixels as vectors ... ie. always get a 1-band image 425// back ... also treat matricies as lists of vectors 426// handy for dE from object difference 427abs_vec x 428 = oo_unary_function abs_vec_op x, is_class x 429 = abs_vec_image x, is_image x 430 = abs_vec_cmplx x, is_complex x 431 = abs_vec_num x, is_real x 432 = abs_vec_list x, is_real_list x 433 = mean (map abs_vec_list x), is_matrix x 434 = error (_ "bad arguments to " ++ "abs_vec") 435{ 436 abs_vec_op = Operator "abs_vec" 437 abs_vec Operator_type.COMPOUND false; 438 439 abs_vec_list l = (sum (map square l)) ** 0.5; 440 441 abs_vec_num n 442 = n, n >= 0 443 = -n; 444 445 abs_vec_cmplx c = ((re c)**2 + (im c)**2) ** 0.5; 446 447 abs_vec_image im 448 = (sum (map square (bandsplit im))) ** 0.5; 449} 450 451transpose x 452 = oo_unary_function transpose_op x, is_class x 453 = transpose_image x, is_image x 454 = transpose_list x, is_listof is_list x 455 = error (_ "bad arguments to " ++ "transpose") 456{ 457 transpose_op = Operator "transpose" 458 transpose Operator_type.COMPOUND_REWRAP false; 459 460 transpose_list l 461 = [], l' == [] 462 = (map hd l') : (transpose_list (map tl l')) 463 { 464 l' = takewhile (not_equal []) l; 465 } 466 467 transpose_image = im_flipver @ im_rot270; 468} 469 470rot45 x 471 = oo_unary_function rot45_op x, is_class x 472 = error "rot45 image: not implemented", is_image x 473 = error (_ "bad arguments to " ++ "rot45") 474{ 475 rot45_op = Operator "rot45" 476 rot45_object Operator_type.COMPOUND_REWRAP false; 477 478 rot45_object x 479 = rot45_matrix x, is_odd_square_matrix x 480 = error "rot45 image: not implemented", is_image x 481 = error (_ "bad arguments to " ++ "rot45"); 482 483 // slow, but what the heck 484 rot45_matrix l = (im_rotate_dmask45 (Matrix l)).value; 485} 486 487// apply an image function to a [[real]] ... matrix is converted to a 1 band 488// image for processing 489apply_matrix_as_image fn m 490 = (get_value @ im_vips2mask @ fn @ im_mask2vips @ Matrix) m; 491 492// a general image/matrix operation where the mat version is most easily done 493// by converting mat->image->mat 494apply_matim_operation name fn x 495 = oo_unary_function class_op x, is_class x 496 = fn x, is_image x 497 = apply_matrix_as_image fn x, is_matrix x 498 = error (_ "bad arguments to " ++ name) 499{ 500 class_op = Operator name 501 (apply_matim_operation name fn) Operator_type.COMPOUND_REWRAP false; 502} 503 504rot90 = apply_matim_operation "rot90" im_rot90; 505rot180 = apply_matim_operation "rot180" im_rot180; 506rot270 = apply_matim_operation "rot270" im_rot270; 507rotquad = apply_matim_operation "rotquad" im_rotquad; 508fliplr = apply_matim_operation "fliplr" im_fliphor; 509fliptb = apply_matim_operation "flipud" im_flipver; 510 511image_set_type type x 512 = oo_unary_function image_set_type_op x, is_class x 513 = im_copy_set x (to_real type) 514 (get_header "Xres" x) (get_header "Yres" x) 515 (get_header "Xoffset" x) (get_header "Yoffset" x), 516 is_image x 517 = error (_ "bad arguments to " ++ "image_set_type:" ++ 518 print type ++ " " ++ print x) 519{ 520 image_set_type_op = Operator "image_set_type" 521 (image_set_type type) Operator_type.COMPOUND_REWRAP false; 522} 523 524image_set_origin xoff yoff x 525 = oo_unary_function image_set_origin_op x, is_class x 526 = im_copy_set x 527 (get_header "Type" x) 528 (get_header "Xres" x) (get_header "Yres" x) 529 (to_real xoff) (to_real yoff), 530 is_image x 531 = error (_ "bad arguments to " ++ "image_set_origin") 532{ 533 image_set_origin_op = Operator "image_set_origin" 534 (image_set_origin xoff yoff) 535 Operator_type.COMPOUND_REWRAP false; 536} 537 538cache tile_width tile_height max_tiles x 539 = oo_unary_function cache_op x, is_class x 540 = im_cache x (to_real tile_width) (to_real tile_height) 541 (to_real max_tiles), is_image x 542 = error (_ "bad arguments to " ++ "cache") 543{ 544 cache_op = Operator "cache" 545 (cache tile_width tile_height max_tiles) 546 Operator_type.COMPOUND_REWRAP false; 547} 548 549tile across down x 550 = oo_unary_function tile_op x, is_class x 551 = im_replicate x (to_real across) (to_real down), is_image x 552 = error (_ "bad arguments to " ++ "tile") 553{ 554 tile_op = Operator "tile" 555 (tile across down) Operator_type.COMPOUND_REWRAP false; 556} 557 558grid tile_height across down x 559 = oo_unary_function grid_op x, is_class x 560 = im_grid x (to_real tile_height) (to_real across) (to_real down), 561 is_image x 562 = error (_ "bad arguments to " ++ "grid") 563{ 564 grid_op = Operator "grid" 565 (grid tile_height across down) Operator_type.COMPOUND_REWRAP false; 566} 567 568max_pair a b 569 = a, a > b 570 = b; 571 572min_pair a b 573 = a, a < b 574 = b; 575 576range min value max = min_pair max (max_pair min value); 577 578max x 579 = oo_unary_function max_op x, is_class x 580 = im_max x, is_image x 581 = max_list x, is_list x 582 = x, is_number x 583 = error (_ "bad arguments to " ++ "max") 584{ 585 max_op = Operator "max" max Operator_type.COMPOUND false; 586 587 max_list x 588 = error "max []", x == [] 589 = foldr1 max_pair x, is_real_list x 590 = foldr1 max_pair (map max_list x), is_list x 591 = max x; 592} 593 594min x 595 = oo_unary_function min_op x, is_class x 596 = im_min x, is_image x 597 = min_list x, is_list x 598 = x, is_number x 599 = error (_ "bad arguments to " ++ "min") 600{ 601 min_op = Operator "min" min Operator_type.COMPOUND false; 602 603 min_list x 604 = error "min []", x == [] 605 = foldr1 min_pair x, is_real_list x 606 = foldr1 min_pair (map min_list x), is_list x 607 = min x; 608} 609 610maxpos x 611 = oo_unary_function maxpos_op x, is_class x 612 = im_maxpos x, is_image x 613 = maxpos_matrix x, is_matrix x 614 = maxpos_list x, is_list x 615 = error (_ "bad arguments to " ++ "maxpos") 616{ 617 maxpos_op = Operator "maxpos" maxpos Operator_type.COMPOUND false; 618 619 maxpos_matrix m 620 = (-1, -1), m == [[]] 621 = (indexes?row, row) 622 { 623 max_value = max (Matrix m); 624 indexes = map (index (equal max_value)) m; 625 row = index (not_equal (-1)) indexes; 626 } 627 628 maxpos_list l 629 = -1, l == [] 630 = index (equal (max l)) l; 631} 632 633minpos x 634 = oo_unary_function minpos_op x, is_class x 635 = im_minpos x, is_image x 636 = minpos_matrix x, is_matrix x 637 = minpos_list x, is_list x 638 = error (_ "bad arguments to " ++ "minpos") 639{ 640 minpos_op = Operator "minpos" minpos Operator_type.COMPOUND false; 641 642 minpos_matrix m 643 = (-1, -1), m == [[]] 644 = (indexes?row, row) 645 { 646 min_value = min (Matrix m); 647 indexes = map (index (equal min_value)) m; 648 row = index (not_equal (-1)) indexes; 649 } 650 651 minpos_list l 652 = -1, l == [] 653 = index (equal (min l)) l; 654} 655 656stats x 657 = oo_unary_function stats_op x, is_class x 658 = im_stats x, is_image x 659 = im_stats (to_image x).value, is_matrix x 660 = error (_ "bad arguments to " ++ "stats") 661{ 662 stats_op = Operator "stats" 663 stats Operator_type.COMPOUND false; 664} 665 666e = 2.7182818284590452354; 667 668pi = 3.14159265358979323846; 669 670rad d = 2 * pi * (d / 360); 671 672deg r = 360 * r / (2 * pi); 673 674sign x 675 = oo_unary_function sign_op x, is_class x 676 = im_sign x, is_image x 677 = sign_cmplx x, is_complex x 678 = sign_num x, is_real x 679 = error (_ "bad arguments to " ++ "sign") 680{ 681 sign_op = Operator "sign" sign Operator_type.COMPOUND_REWRAP false; 682 683 sign_num n 684 = 0, n == 0 685 = 1, n > 0 686 = -1; 687 688 sign_cmplx c 689 = (0, 0), mod == 0 690 = (re c / mod, im c / mod) 691 { 692 mod = abs c; 693 } 694} 695 696rint x 697 = oo_unary_function rint_op x, is_class x 698 = im_rint x, is_image x 699 = rint_value x, is_number x 700 = error (_ "bad arguments to " ++ "rint") 701{ 702 rint_op = Operator "rint" rint Operator_type.ARITHMETIC false; 703 704 rint_value x 705 = (int) (x + 0.5), x > 0 706 = (int) (x - 0.5); 707} 708 709scale x 710 = oo_unary_function scale_op x, is_class x 711 = (unsigned char) x, is_number x 712 = im_scale x, is_image x 713 = scale_list x, is_real_list x || is_matrix x 714 = error (_ "bad arguments to " ++ "scale") 715{ 716 scale_op = Operator "scale" scale Operator_type.COMPOUND_REWRAP false; 717 718 scale_list l 719 = apply_scale s o l 720 { 721 mn = find_limit min_pair l; 722 mx = find_limit max_pair l; 723 s = 255.0 / (mx - mn); 724 o = -(mn * s); 725 } 726 727 find_limit fn l 728 = find_limit fn (map (find_limit fn) l), is_listof is_list l 729 = foldr1 fn l; 730 731 apply_scale s o x 732 = x * s + o, is_number x 733 = map (apply_scale s o) x; 734} 735 736scaleps x 737 = oo_unary_function scale_op x, is_class x 738 = im_scaleps x, is_image x 739 = error (_ "bad arguments to " ++ "scale") 740{ 741 scale_op = Operator "scaleps" 742 scaleps Operator_type.COMPOUND_REWRAP false; 743} 744 745fwfft x 746 = oo_unary_function fwfft_op x, is_class x 747 = im_fwfft x, is_image x 748 = error (_ "bad arguments to " ++ "fwfft") 749{ 750 fwfft_op = Operator "fwfft" 751 fwfft Operator_type.COMPOUND_REWRAP false; 752} 753 754invfft x 755 = oo_unary_function invfft_op x, is_class x 756 = im_invfftr x, is_image x 757 = error (_ "bad arguments to " ++ "invfft") 758{ 759 invfft_op = Operator "invfft" 760 invfft Operator_type.COMPOUND_REWRAP false; 761} 762 763falsecolour x 764 = oo_unary_function falsecolour_op x, is_class x 765 = image_set_type Image_type.sRGB (im_falsecolour x), is_image x 766 = error (_ "bad arguments to " ++ "falsecolour") 767{ 768 falsecolour_op = Operator "falsecolour" 769 falsecolour Operator_type.COMPOUND_REWRAP false; 770} 771 772polar x 773 = oo_unary_function polar_op x, is_class x 774 = im_c2amph x, is_image x 775 = polar_cmplx x, is_complex x 776 = error (_ "bad arguments to " ++ "polar") 777{ 778 polar_op = Operator "polar" polar Operator_type.COMPOUND false; 779 780 polar_cmplx r 781 = (l, a) 782 { 783 a 784 = 270, x == 0 && y < 0 785 = 90, x == 0 && y >= 0 786 = 360 + atan (y / x), x > 0 && y < 0 787 = atan (y / x), x > 0 && y >= 0 788 = 180 + atan (y / x); 789 790 l = (x ** 2 + y ** 2) ** 0.5; 791 792 x = re r; 793 y = im r; 794 } 795} 796 797rectangular x 798 = oo_unary_function rectangular_op x, is_class x 799 = im_c2rect x, is_image x 800 = rectangular_cmplx x, is_complex x 801 = error (_ "bad arguments to " ++ "rectangular") 802{ 803 rectangular_op = Operator "rectangular" 804 rectangular Operator_type.COMPOUND false; 805 806 rectangular_cmplx p 807 = (x, y) 808 { 809 l = re p; 810 a = im p; 811 812 x = l * cos a; 813 y = l * sin a; 814 } 815} 816 817// we can't use colour_unary: that likes 3 band only 818recomb matrix x 819 = oo_unary_function recomb_op x, is_class x 820 = im_recomb x matrix, is_image x 821 = recomb_real_list x, is_real_list x 822 = map recomb_real_list x, is_matrix x 823 = error (_ "bad arguments to " ++ "recomb") 824{ 825 // COMPOUND_REWRAP ... signal to the colour class to go to image and 826 // back 827 recomb_op = Operator "recomb" 828 (recomb matrix) Operator_type.COMPOUND_REWRAP false; 829 830 // process [1,2,3 ..] as an image 831 recomb_real_list l 832 = (to_matrix im').value?0 833 { 834 im = (float) (to_image (Vector l)).value; 835 im' = recomb matrix im; 836 } 837} 838 839extract_area x y w h obj 840 = oo_unary_function extract_area_op obj, is_class obj 841 = im_extract_area obj x' y' w' h', is_image obj 842 = map (extract_range x' w') (extract_range y' h' obj), is_matrix obj 843 = error (_ "bad arguments to " ++ "extract_area") 844{ 845 x' = to_real x; 846 y' = to_real y; 847 w' = to_real w; 848 h' = to_real h; 849 850 extract_area_op = Operator "extract_area" (extract_area x y w h) 851 Operator_type.COMPOUND_REWRAP false; 852 853 extract_range from length list 854 = (take length @ drop from) list; 855} 856 857extract_band b obj = subscript obj b; 858 859extract_row y obj 860 = oo_unary_function extract_row_op obj, is_class obj 861 = extract_area 0 y' (get_width obj) 1 obj, is_image obj 862 = [obj?y'], is_matrix obj 863 = error (_ "bad arguments to " ++ "extract_row") 864{ 865 y' = to_real y; 866 867 extract_row_op = Operator "extract_row" (extract_row y) 868 Operator_type.COMPOUND_REWRAP false; 869} 870 871extract_column x obj 872 = oo_unary_function extract_column_op obj, is_class obj 873 = extract_area x' 0 1 height obj, is_image obj 874 = map (\row [row?x']) obj, is_matrix obj 875 = error (_ "bad arguments to " ++ "extract_column") 876{ 877 x' = to_real x; 878 height = get_header "Ysize" obj; 879 880 extract_column_op = Operator "extract_column" (extract_column x) 881 Operator_type.COMPOUND_REWRAP false; 882} 883 884blend cond in1 in2 885 = oo_binary_function blend_op cond [in1,in2], is_class cond 886 = im_blend (get_image cond) (get_image in1) (get_image in2), 887 has_image cond && has_image in1 && has_image in2 888 = error (_ "bad arguments to " ++ "blend") 889{ 890 blend_op = Operator "blend" 891 blend_obj Operator_type.COMPOUND_REWRAP false; 892 893 blend_obj cond x 894 = blend_result_image 895 { 896 [then_part, else_part] = x; 897 898 // get things about our output from inputs in this order 899 objects = [then_part, else_part, cond]; 900 901 // properties of our output image 902 target_width = get_member_list has_width get_width objects; 903 target_height = get_member_list has_height get_height objects; 904 target_bands = get_member_list has_bands get_bands objects; 905 target_format = get_member_list has_format get_format objects; 906 target_type = get_member_list has_type get_type objects; 907 908 to_image x 909 = x, is_image x 910 = x.value, is_Image x 911 = black + x 912 { 913 black = im_black target_width target_height target_bands; 914 } 915 916 [then_image, else_image] = 917 map (clip2fmt target_format @ to_image) [then_part, else_part]; 918 [c, t, e] = size_alike [cond, then_image, else_image]; 919 920 blend_result_image = image_set_type target_type (im_blend c t e); 921 } 922} 923 924insert x y small big 925 = oo_binary_function insert_op small big, is_class small 926 = oo_binary'_function insert_op small big, is_class big 927 = im_insert big' small' (to_real x) (to_real y), 928 is_image small && is_image big 929 = error (_ "bad arguments to " ++ "insert") 930{ 931 insert_op = Operator "insert" 932 (insert x y) Operator_type.COMPOUND_REWRAP false; 933 [small', big'] = (formats_alike @ bands_alike) [small, big]; 934} 935 936insert_noexpand x y small big 937 = oo_binary_function insert_noexpand_op small big, is_class small 938 = oo_binary'_function insert_noexpand_op small big, is_class big 939 = im_insert_noexpand big' small' (to_real x) (to_real y), 940 is_image small && is_image big 941 = error (_ "bad arguments to " ++ "insert_noexpand") 942{ 943 insert_noexpand_op = Operator "insert_noexpand" 944 (insert_noexpand x y) Operator_type.COMPOUND_REWRAP false; 945 [small', big'] = (formats_alike @ bands_alike) [small, big]; 946} 947 948measure x y w h u v image 949 = oo_unary_function measure_op image, is_class image 950 = im_measure image 951 (to_real x) (to_real y) (to_real w) (to_real h) 952 (to_real u) (to_real v), 953 is_image image 954 = error (_ "bad arguments to " ++ "measure") 955{ 956 measure_op = Operator "measure" 957 (measure x y w h u v) Operator_type.COMPOUND_REWRAP false; 958} 959 960extract_bands b n obj 961 = oo_unary_function extract_bands_op obj, is_class obj 962 = im_extract_bands obj (to_real b) (to_real n), is_image obj 963 = error (_ "bad arguments to " ++ "extract_bands") 964{ 965 extract_bands_op = Operator "extract_bands" 966 (extract_bands b n) Operator_type.COMPOUND_REWRAP false; 967} 968 969invert x 970 = oo_unary_function invert_op x, is_class x 971 = im_invert x, is_image x 972 = 255 - x, is_real x 973 = error (_ "bad arguments to " ++ "invert") 974{ 975 invert_op = Operator "invert" invert Operator_type.COMPOUND false; 976} 977 978transform ipol wrap params image 979 = oo_unary_function transform_op image, is_class image 980 = im_transform image 981 (to_matrix params) (to_real ipol) (to_real wrap), is_image image 982 = error (_ "bad arguments to " ++ "transform") 983{ 984 transform_op = Operator "transform" 985 (transform ipol wrap params) 986 Operator_type.COMPOUND_REWRAP false; 987} 988 989transform_search max_error max_iterations order ipol wrap sample reference 990 = oo_binary_function transform_search_op sample reference, is_class sample 991 = oo_binary'_function transform_search_op sample reference, 992 is_class reference 993 = im_transform_search sample reference 994 (to_real max_error) (to_real max_iterations) (to_real order) 995 (to_real ipol) (to_real wrap), 996 is_image sample && is_image reference 997 = error (_ "bad arguments to " ++ "transform_search") 998{ 999 transform_search_op = Operator "transform_search" 1000 (transform_search max_error max_iterations order ipol wrap) 1001 Operator_type.COMPOUND false; 1002} 1003 1004rotate angle image 1005 = oo_binary_function rotate_op angle image, is_class angle 1006 = oo_binary'_function rotate_op angle image, is_class image 1007 = im_similarity image (cos angle) (sin angle) 0 0, 1008 is_real angle && is_image image 1009 = error (_ "bad arguments to " ++ "rotate") 1010{ 1011 rotate_op = Operator "rotate" 1012 rotate Operator_type.COMPOUND_REWRAP false; 1013} 1014 1015matrix_binary fn a b 1016 = itom (fn (mtoi a) (mtoi b)) 1017{ 1018 mtoi x = im_mask2vips (Matrix x); 1019 itom x = (im_vips2mask x).value; 1020} 1021 1022join_lr a b 1023 = oo_binary_function join_lr_op a b, is_class a 1024 = oo_binary'_function join_lr_op a b, is_class b 1025 = join a b 1026{ 1027 join_lr_op = Operator "join_lr" 1028 join Operator_type.COMPOUND_REWRAP false; 1029 1030 join a b 1031 = join_im a b, is_image a && is_image b 1032 = matrix_binary join_im a b, 1033 is_matrix a && is_matrix b 1034 = error (_ "bad arguments to " ++ "join_lr"); 1035 1036 join_im a b = insert (get_width a) 0 b a; 1037} 1038 1039join_tb a b 1040 = oo_binary_function join_tb_op a b, is_class a 1041 = oo_binary'_function join_tb_op a b, is_class b 1042 = join a b 1043{ 1044 join_tb_op = Operator "join_tb" 1045 join Operator_type.COMPOUND_REWRAP false; 1046 1047 join a b 1048 = join_im a b, is_image a && is_image b 1049 = matrix_binary join_im a b, 1050 is_matrix a && is_matrix b 1051 = error (_ "bad arguments to " ++ "join_tb"); 1052 1053 join_im a b = insert 0 (get_height a) b a; 1054} 1055 1056conj x 1057 = oo_unary_function conj_op x, is_class x 1058 = (re x, -im x), 1059 is_complex x || 1060 (is_image x && format == Image_format.COMPLEX) || 1061 (is_image x && format == Image_format.DPCOMPLEX) 1062 // assume it's some sort of real 1063 = x 1064{ 1065 format = get_header "BandFmt" x; 1066 conj_op = Operator "conj" conj Operator_type.COMPOUND false; 1067} 1068 1069clip2fmt format image 1070 = oo_unary_function clip2fmt_op image, is_class image 1071 = im_clip2fmt image (to_real format), is_image image 1072 = error (_ "bad arguments to " ++ "clip2fmt") 1073{ 1074 clip2fmt_op = Operator "clip2fmt" 1075 (clip2fmt format) Operator_type.COMPOUND_REWRAP false; 1076} 1077 1078embed type x y w h im 1079 = oo_unary_function embed_op im, is_class im 1080 = im_embed im (to_real type) 1081 (to_real x) (to_real y) (to_real w) (to_real h), is_image im 1082 = error (_ "bad arguments to " ++ "embed") 1083{ 1084 embed_op = Operator "embed" 1085 (embed type x y w h) Operator_type.COMPOUND_REWRAP false; 1086} 1087 1088/* Morph a mask with a [[real]] matrix ... turn m2 into an image, morph it 1089 * with m1, turn it back to a matrix again. 1090 */ 1091_morph_2_masks fn m1 m2 1092 = m'' 1093{ 1094 image = (unsigned char) im_mask2vips (Matrix m2); 1095 m2_width = get_width image; 1096 m2_height = get_height image; 1097 1098 // need to embed m2 in an image large enough for us to be able to 1099 // position m1 all around the edges, with a 1 pixel overlap 1100 image' = embed 0 1101 (m1.width / 2) (m1.height / 2) 1102 (m2_width + (m1.width - 1)) (m2_height + (m1.height - 1)) 1103 image; 1104 1105 // morph! 1106 image'' = fn m1 image'; 1107 1108 // back to mask 1109 m' = im_vips2mask ((double) image''); 1110 1111 // Turn 0 in output to 128 (don't care). 1112 m'' 1113 = map (map fn) m'.value 1114 { 1115 fn a 1116 = 128, a == 0; 1117 = a; 1118 } 1119} 1120 1121dilate mask image 1122 = oo_unary_function dilate_op image, is_class image 1123 = im_dilate image (to_matrix mask), is_image image 1124 = error (_ "bad arguments to " ++ "dilate") 1125{ 1126 dilate_op = Operator "dilate" 1127 dilate_object Operator_type.COMPOUND_REWRAP false; 1128 1129 dilate_object x 1130 = _morph_2_masks dilate mask x, is_matrix x 1131 = dilate mask x; 1132} 1133 1134erode mask image 1135 = oo_unary_function erode_op image, is_class image 1136 = im_erode image (to_matrix mask), is_image image 1137 = error (_ "bad arguments to " ++ "erode") 1138{ 1139 erode_op = Operator "erode" 1140 erode_object Operator_type.COMPOUND_REWRAP false; 1141 1142 erode_object x 1143 = _morph_2_masks erode mask x, is_matrix x 1144 = erode mask x; 1145} 1146 1147conv mask image 1148 = oo_unary_function conv_op image, is_class image 1149 = im_conv image (to_matrix mask), is_image image 1150 = error (_ "bad arguments to " ++ "conv") 1151{ 1152 conv_op = Operator "conv" 1153 (conv mask) Operator_type.COMPOUND_REWRAP false; 1154} 1155 1156convsep mask image 1157 = oo_unary_function convsep_op image, is_class image 1158 = im_convsep image (to_matrix mask), is_image image 1159 = error (_ "bad arguments to " ++ "convsep") 1160{ 1161 convsep_op = Operator "convsep" 1162 (convsep mask) Operator_type.COMPOUND_REWRAP false; 1163} 1164 1165convsepf mask image 1166 = oo_unary_function convsepf_op image, is_class image 1167 = im_convsepf image (to_matrix mask), is_image image 1168 = error (_ "bad arguments to " ++ "convsepf") 1169{ 1170 convsepf_op = Operator "convsepf" 1171 (convsepf mask) Operator_type.COMPOUND_REWRAP false; 1172} 1173 1174rank w h n image 1175 = oo_unary_function rank_op image, is_class image 1176 = im_rank image (to_real w) (to_real h) (to_real n), is_image image 1177 = error (_ "bad arguments to " ++ "rank") 1178{ 1179 rank_op = Operator "rank" 1180 (rank w h n) Operator_type.COMPOUND_REWRAP false; 1181} 1182 1183rank_image n x 1184 = rlist x.value, is_Group x 1185 = rlist x, is_list x 1186 = error (_ "bad arguments to " ++ "rank_image") 1187{ 1188 rlist l 1189 = wrapper ranked, has_wrapper 1190 = ranked 1191 { 1192 has_wrapper = has_member_list (has_member "Image") l; 1193 wrapper = get_member_list (has_member "Image") (get_member "Image") l; 1194 ranked = im_rank_image (map get_image l) (to_real n); 1195 } 1196} 1197 1198greyc iterations amplitude sharpness anisotropy alpha 1199 sigma dl da gauss_prec interpolation fast_approx x 1200 = oo_unary_function greyc_op x, is_class x 1201 = greyc_im x, is_image x 1202 = error (_ "bad argument" ++ " (" ++ print x ++ ") to " ++ "greyc") 1203{ 1204 greyc_op = Operator "greyc" (greyc 1205 iterations amplitude sharpness anisotropy alpha 1206 sigma dl da gauss_prec interpolation fast_approx) 1207 Operator_type.COMPOUND_REWRAP false; 1208 greyc_im x = im_greyc x 1209 iterations amplitude sharpness anisotropy alpha 1210 sigma dl da gauss_prec interpolation fast_approx; 1211} 1212 1213greyc_mask iterations amplitude sharpness anisotropy alpha 1214 sigma dl da gauss_prec interpolation fast_approx mask x 1215 = oo_binary_function greyc_mask_op mask x, is_class mask 1216 = oo_binary'_function greyc_mask_op mask x, is_class x 1217 = greyc_im mask x, is_image mask && is_image x 1218 = error (_ "bad arguments" ++ 1219 " (" ++ print mask ++ ", " ++ print x ++ ") " ++ 1220 "to " ++ "greyc") 1221{ 1222 greyc_mask_op = Operator "greyc_mask" (greyc_mask 1223 iterations amplitude sharpness anisotropy alpha 1224 sigma dl da gauss_prec interpolation fast_approx) 1225 Operator_type.COMPOUND_REWRAP false; 1226 greyc_im mask x = im_greyc_mask x mask 1227 iterations amplitude sharpness anisotropy alpha 1228 sigma dl da gauss_prec interpolation fast_approx; 1229} 1230 1231// find the correlation surface for a small image within a big one 1232correlate small big 1233 = oo_binary_function correlate_op small big, is_class small 1234 = oo_binary'_function correlate_op small big, is_class big 1235 = im_spcor big small, is_image small && is_image big 1236 = error (_ "bad arguments to " ++ "correlate") 1237{ 1238 correlate_op = Operator "correlate" 1239 correlate Operator_type.COMPOUND_REWRAP false; 1240} 1241 1242// just sum-of-squares-of-differences 1243correlate_fast small big 1244 = oo_binary_function correlate_fast_op small big, is_class small 1245 = oo_binary'_function correlate_fast_op small big, is_class big 1246 = im_fastcor big small, is_image small && is_image big 1247 = error (_ "bad arguments to " ++ "correlate_fast") 1248{ 1249 correlate_fast_op = Operator "correlate_fast" 1250 correlate_fast Operator_type.COMPOUND_REWRAP false; 1251} 1252 1253// set Type, wrap as Plot_hist if it's a class 1254hist_tag x 1255 = oo_unary_function hist_tag_op x, is_class x 1256 = image_set_type Image_type.HISTOGRAM x, is_image x 1257 = error (_ "bad arguments to " ++ "hist_tag") 1258{ 1259 hist_tag_op = Operator "hist_tag" 1260 (Plot_histogram @ hist_tag) Operator_type.COMPOUND false; 1261} 1262 1263hist_find x 1264 = oo_unary_function hist_find_op x, is_class x 1265 = im_histgr x (-1), is_image x 1266 = error (_ "bad arguments to " ++ "hist_find") 1267{ 1268 hist_find_op = Operator "hist_find" 1269 (Plot_histogram @ hist_find) Operator_type.COMPOUND false; 1270} 1271 1272hist_find_nD bins image 1273 = oo_unary_function hist_find_nD_op image, is_class image 1274 = im_histnD image (to_real bins), is_image image 1275 = error (_ "bad arguments to " ++ "hist_find_nD") 1276{ 1277 hist_find_nD_op = Operator "hist_find_nD" 1278 (hist_find_nD bins) Operator_type.COMPOUND_REWRAP false; 1279} 1280 1281hist_map hist image 1282 = oo_binary_function hist_map_op hist image, is_class hist 1283 = oo_binary'_function hist_map_op hist image, is_class image 1284 = im_maplut image hist, is_image hist && is_image image 1285 = error (_ "bad arguments to " ++ "hist_map") 1286{ 1287 // can't use rewrap, as we want to always wrap as image 1288 hist_map_op = Operator "hist_map" 1289 (compose (compose Image) hist_map) Operator_type.COMPOUND false; 1290} 1291 1292hist_cum hist 1293 = oo_unary_function hist_cum_op hist, is_class hist 1294 = im_histcum hist, is_image hist 1295 = error (_ "bad arguments to " ++ "hist_cum") 1296{ 1297 hist_cum_op = Operator "hist_cum" 1298 hist_cum Operator_type.COMPOUND_REWRAP false; 1299} 1300 1301hist_diff hist 1302 = oo_unary_function hist_diff_op hist, is_class hist 1303 = im_histdiff hist, is_image hist 1304 = error (_ "bad arguments to " ++ "hist_diff") 1305{ 1306 hist_diff_op = Operator "hist_diff" 1307 hist_diff Operator_type.COMPOUND_REWRAP false; 1308 1309 im_histdiff h 1310 = (conv (Matrix [[-1, 1]]) @ clip2fmt (fmt (get_format h))) h 1311 { 1312 // up the format so it can represent the whole range of 1313 // possible values from this mask 1314 fmt x 1315 = Image_format.SHORT, 1316 x == Image_format.UCHAR || x == Image_format.CHAR 1317 = Image_format.INT, 1318 x == Image_format.USHORT || x == Image_format.SHORT || 1319 x == Image_format.UINT 1320 = x; 1321 } 1322} 1323 1324hist_norm hist 1325 = oo_unary_function hist_norm_op hist, is_class hist 1326 = im_histnorm hist, is_image hist 1327 = error (_ "bad arguments to " ++ "hist_norm") 1328{ 1329 hist_norm_op = Operator "hist_norm" 1330 hist_norm Operator_type.COMPOUND_REWRAP false; 1331} 1332 1333hist_match in ref 1334 = oo_binary_function hist_match_op in ref, is_class in 1335 = oo_binary'_function hist_match_op in ref, is_class ref 1336 = im_histspec in ref, is_image in && is_image ref 1337 = error (_ "bad arguments to " ++ "hist_match") 1338{ 1339 hist_match_op = Operator "hist_match" 1340 hist_match Operator_type.COMPOUND_REWRAP false; 1341} 1342 1343hist_equalize x = hist_map ((hist_norm @ hist_cum @ hist_find) x) x; 1344 1345hist_equalize_local w h image 1346 = oo_unary_function hist_equalize_local_op image, is_class image 1347 = lhisteq image, is_image image 1348 = error (_ "bad arguments to " ++ "hist_equalize_local") 1349{ 1350 hist_equalize_local_op = Operator "hist_equalize_local" 1351 (hist_equalize_local w h) Operator_type.COMPOUND_REWRAP false; 1352 1353 // loop over bands, if necessary 1354 lhisteq im 1355 = im_lhisteq im (to_real w) (to_real h), get_bands im == 1 1356 = (foldl1 join @ map lhisteq @ bandsplit) im; 1357} 1358 1359// find the threshold below which are percent of the image (percent in [0,1]) 1360// eg. hist_thresh 0.1 x == 12, then x < 12 will light up 10% of the pixels 1361hist_thresh percent image 1362 = x 1363{ 1364 // our own normaliser ... we don't want to norm channels separately 1365 // norm to [0,1] 1366 my_hist_norm h = h / max h; 1367 1368 // normalised cumulative hist 1369 // we sum the channels before we normalise, because we want to treat them 1370 // all the same 1371 h = (my_hist_norm @ sum @ bandsplit @ hist_cum @ hist_find) 1372 image.value; 1373 1374 // threshold that, then use im_profile to search for the x position in the 1375 // histogram 1376 x = mean (im_profile (h > percent) 1); 1377} 1378 1379/* Sometimes useful, despite plotting now being built in, for making 1380 * diagnostic images. 1381 */ 1382hist_plot hist 1383 = oo_unary_function hist_plot_op hist, is_class hist 1384 = im_histplot hist, is_image hist 1385 = error (_ "bad arguments to " ++ "hist_plot") 1386{ 1387 hist_plot_op = Operator "hist_plot" 1388 (Image @ hist_plot) Operator_type.COMPOUND false; 1389} 1390 1391zerox d x 1392 = oo_unary_function zerox_op x, is_class x 1393 = im_zerox x (to_real d), is_image x 1394 = error (_ "bad arguments to " ++ "zerox") 1395{ 1396 zerox_op = Operator "zerox" 1397 (zerox d) Operator_type.COMPOUND_REWRAP false; 1398} 1399 1400resize xfac yfac interp image 1401 = oo_unary_function resize_op image, is_class image 1402 = resize_im image, is_image image 1403 = error (_ "bad arguments to " ++ "resize") 1404{ 1405 resize_op = Operator "resize" 1406 resize_im Operator_type.COMPOUND_REWRAP false; 1407 1408 xfac' = to_real xfac; 1409 yfac' = to_real yfac; 1410 1411 rxfac' = 1 / xfac'; 1412 ryfac' = 1 / yfac'; 1413 1414 resize_im im 1415 // upscale by integer factor, nearest neighbour 1416 = im_zoom im xfac' yfac', 1417 is_int xfac' && is_int yfac' && 1418 xfac' >= 1 && yfac' >= 1 && 1419 interp == Interpolate.NEAREST_NEIGHBOUR 1420 1421 // downscale by integer factor, nearest neighbour 1422 = im_subsample im rxfac' ryfac', 1423 is_int rxfac' && is_int ryfac' && 1424 rxfac' >= 1 && ryfac' >= 1 && 1425 interp == Interpolate.NEAREST_NEIGHBOUR 1426 1427 // upscale by any factor, nearest neighbour 1428 // can't really do this right ... upscale by integer part, then 1429 // bilinear to exact size 1430 = scale xg?1 yg?1 (im_zoom im xg?0 yg?0), 1431 xfac' >= 1 && yfac' >= 1 && 1432 interp == Interpolate.NEAREST_NEIGHBOUR 1433 1434 // downscale by any factor, nearest neighbour 1435 // can't really do this right ... downscale by integer part, 1436 // then bilinear to exact size 1437 = scale xs?1 ys?1 (im_subsample im xs?0 ys?0), 1438 rxfac' >= 1 && ryfac' >= 1 && 1439 interp == Interpolate.NEAREST_NEIGHBOUR 1440 1441 // upscale by any factor, bilinear 1442 = scale xfac' yfac' im, 1443 xfac' >= 1 && yfac' >= 1 && 1444 interp == Interpolate.BILINEAR 1445 1446 // downscale by any factor, bilinear 1447 // block shrink by integer factor, then bilinear resample to 1448 // exact 1449 = scale xs?1 ys?1 (im_shrink im xs?0 ys?0), 1450 rxfac' >= 1 && ryfac' >= 1 && 1451 interp == Interpolate.BILINEAR 1452 1453 = error ("resize: unimplemented argument combination:\n" ++ 1454 " xfac = " ++ print xfac' ++ "\n" ++ 1455 " yfac = " ++ print yfac' ++ "\n" ++ 1456 " interp = " ++ print interp ++ " (" ++ 1457 Interpolate.names.lookup 1 0 interp ++ ")") 1458 { 1459 // convert a float scale to integer plus fraction 1460 // eg. scale by 2.5 becomes [2, 1.25] (x * 2.5 == x * 2 * 1.25) 1461 break f = [floor f, f / floor f]; 1462 1463 // same, but for downsizing ... turn a float scale which is less than 1464 // 1 into an int shrink and a float scale 1465 1466 // complicated: the int shrink may round the size down (eg. imagine 1467 // subsampling a 11 pixel wide image by 3, we'd get a 3 pixel wide 1468 // image, not a 3.666 pixel wide image), so pass in the size of image 1469 // we are operating on and adjust for any rounding 1470 1471 // so ... x is (eg.) 467, f is (eg. 128/467, about 0.274) 1472 rbreak x f 1473 = [int_shrink, float_resample] 1474 { 1475 // the size of image we are aiming for after the combined int and 1476 // float resample 1477 x' = x * f; 1478 1479 // int part 1480 int_shrink = floor (1 / f); 1481 1482 // size after int shrink 1483 x'' = floor (x / int_shrink); 1484 1485 // therefore what we need for the float part 1486 float_resample = x' / x''; 1487 } 1488 1489 width = get_width im; 1490 height = get_height im; 1491 1492 // grow and shrink factors 1493 xg = break xfac'; 1494 yg = break yfac'; 1495 xs = rbreak width xfac'; 1496 ys = rbreak height yfac'; 1497 1498 // binlinear resize 1499 scale xfac yfac im 1500 = im_affine im 1501 xfac 0 0 yfac 1502 0 0 1503 0 0 1504 (rint (get_width im * xfac)) 1505 (rint (get_height im * yfac)); 1506 } 1507} 1508 1509sharpen radius x1 y2 y3 m1 m2 in 1510 = oo_unary_function sharpen_op in, is_class in 1511 = im_sharpen in (to_real radius) 1512 (to_real x1) (to_real y2) (to_real y3) 1513 (to_real m1) (to_real m2), is_image in 1514 = error (_ "bad arguments to " ++ "sharpen") 1515{ 1516 sharpen_op = Operator "sharpen" 1517 (sharpen radius x1 y2 y3 m1 m2) 1518 Operator_type.COMPOUND_REWRAP false; 1519} 1520 1521tone_analyse s m h sa ma ha in 1522 = oo_unary_function tone_analyse_op in, is_class in 1523 = im_tone_analyse in 1524 (to_real s) (to_real m) (to_real h) 1525 (to_real sa) (to_real ma) (to_real ha), is_image in 1526 = error (_ "bad arguments to " ++ "tone_analyse") 1527{ 1528 tone_analyse_op = Operator "tone_analyse" 1529 (Plot_histogram @ tone_analyse s m h sa ma ha) 1530 Operator_type.COMPOUND false; 1531} 1532 1533tone_map hist image 1534 = oo_binary_function tone_map_op hist image, is_class hist 1535 = oo_binary'_function tone_map_op hist image, is_class image 1536 = im_tone_map image hist, is_image hist && is_image image 1537 = error (_ "bad arguments to " ++ "tone_map") 1538{ 1539 tone_map_op = Operator "tone_map" 1540 tone_map Operator_type.COMPOUND_REWRAP false; 1541} 1542 1543tone_build fmt b w s m h sa ma ha 1544 = (Plot_histogram @ clip2fmt fmt) 1545 (im_tone_build_range mx mx 1546 (to_real b) (to_real w) 1547 (to_real s) (to_real m) (to_real h) 1548 (to_real sa) (to_real ma) (to_real ha)) 1549{ 1550 mx = Image_format.maxval fmt; 1551} 1552 1553icc_export depth profile intent in 1554 = oo_unary_function icc_export_op in, is_class in 1555 = im_icc_export_depth in 1556 (to_real depth) (expand profile) (to_real intent), is_image in 1557 = error (_ "bad arguments to " ++ "icc_export") 1558{ 1559 icc_export_op = Operator "icc_export" 1560 (icc_export depth profile intent) 1561 Operator_type.COMPOUND_REWRAP false; 1562} 1563 1564icc_import profile intent in 1565 = oo_unary_function icc_import_op in, is_class in 1566 = im_icc_import in 1567 (expand profile) (to_real intent), is_image in 1568 = error (_ "bad arguments to " ++ "icc_import") 1569{ 1570 icc_import_op = Operator "icc_import" 1571 (icc_import profile intent) 1572 Operator_type.COMPOUND_REWRAP false; 1573} 1574 1575icc_import_embedded intent in 1576 = oo_unary_function icc_import_embedded_op in, is_class in 1577 = im_icc_import_embedded in (to_real intent), is_image in 1578 = error (_ "bad arguments to " ++ "icc_import_embedded") 1579{ 1580 icc_import_embedded_op = Operator "icc_import_embedded" 1581 (icc_import_embedded intent) 1582 Operator_type.COMPOUND_REWRAP false; 1583} 1584 1585icc_transform in_profile out_profile intent in 1586 = oo_unary_function icc_transform_op in, is_class in 1587 = im_icc_transform in 1588 (expand in_profile) (expand out_profile) 1589 (to_real intent), is_image in 1590 = error (_ "bad arguments to " ++ "icc_transform") 1591{ 1592 icc_transform_op = Operator "icc_transform" 1593 (icc_transform in_profile out_profile intent) 1594 Operator_type.COMPOUND_REWRAP false; 1595} 1596 1597icc_ac2rc profile in 1598 = oo_unary_function icc_ac2rc_op in, is_class in 1599 = im_icc_ac2rc in (expand profile), is_image in 1600 = error (_ "bad arguments to " ++ "icc_ac2rc") 1601{ 1602 icc_ac2rc_op = Operator "icc_ac2rc" 1603 (icc_ac2rc profile) 1604 Operator_type.COMPOUND_REWRAP false; 1605} 1606 1607flood_blob x y v image 1608 = oo_unary_function flood_blob_op image, is_class image 1609 = im_flood_blob_copy image (to_real x) (to_real y) v, is_image image 1610 = error (_ "bad arguments to " ++ "flood_blob") 1611{ 1612 flood_blob_op = Operator "flood_blob" 1613 (flood_blob x y v) Operator_type.COMPOUND_REWRAP false; 1614} 1615 1616print_base base in 1617 = oo_unary_function print_base_op in, is_class in 1618 = map (print_base base) in, is_list in 1619 = print_base_real, is_real in 1620 = error (_ "bad arguments to " ++ "print_base") 1621{ 1622 print_base_op 1623 = Operator "print_base" (print_base base) Operator_type.COMPOUND false; 1624 1625 print_base_real 1626 = error "print_base: bad base", base < 2 || base > 16 1627 = "0", in < 0 || chars == [] 1628 = reverse chars 1629 { 1630 digits = map (\x x % base) 1631 (takewhile (not_equal 0) 1632 (iterate (\x idivide x base) in)); 1633 chars = map tohd digits; 1634 1635 tohd x 1636 = (char) ((int) '0' + x), x < 10 1637 = (char) ((int) 'A' + (x - 10)); 1638 } 1639} 1640 1641/* id x: the identity function 1642 * 1643 * id :: * -> * 1644 */ 1645id x = x; 1646 1647/* const x y: junk y, return x 1648 * 1649 * (const 3) is the function that always returns 3. 1650 * const :: * -> ** -> * 1651 */ 1652const x y = x; 1653 1654/* converse fn a b: swap order of args to fn 1655 * 1656 * converse fn a b == fn b a 1657 * converse :: (* -> ** -> ***) -> ** -> * -> *** 1658 */ 1659converse fn a b = fn b a; 1660 1661/* fix fn x: find the fixed point of a function 1662 */ 1663fix fn x = limit (iterate fn x); 1664 1665/* until pred fn n: apply fn to n until pred succeeds; return that value 1666 * 1667 * until (more 1000) (multiply 2) 1 = 1024 1668 * until :: (* -> bool) -> (* -> *) -> * -> * 1669 */ 1670until pred fn n 1671 = n, pred n 1672 = until pred fn (fn n); 1673 1674/* Infinite list of primes. 1675 */ 1676primes 1677 = 1 : (sieve [2 ..]) 1678{ 1679 sieve l = hd l : sieve (filter (nmultiple (hd l)) (tl l)); 1680 nmultiple n x = x / n != (int) (x / n); 1681} 1682 1683/* Map an n-ary function (pass the args as a list) over groups of objects. 1684 * The objects can be single objects or groups. If more than one 1685 * object is a group, we iterate for the length of the smallest group. 1686 * Don't loop over plain lists, since we want (eg.) "mean [1, 2, 3]" to work. 1687 * Treat [] as no-value --- ie. if any of the inputs are [] we put [] into the 1688 * output and don't apply the function. 1689 1690 copy-pasted into _types, keep in sync 1691 1692 */ 1693map_nary fn args 1694 = fn args, groups == [] 1695 = Group (map process [0, 1 .. shortest - 1]) 1696{ 1697 // find all the group arguments 1698 groups = filter is_Group args; 1699 1700 // what's the length of the shortest group arg? 1701 shortest = foldr1 min_pair (map (len @ get_value) groups); 1702 1703 // process index n ... pull that member from each argument 1704 // recurse to handle application, so we work for nested groups too 1705 process n 1706 = NULL, any (map (is_noval n) args) 1707 = map_nary fn (map (extract n) args) 1708 { 1709 extract n arg 1710 = arg.value?n, is_Group arg 1711 = arg; 1712 1713 is_noval n arg = is_Group arg && arg.value?n == NULL; 1714 } 1715} 1716 1717/* Map a 1-ary function over an object. 1718 */ 1719map_unary fn a = map_nary (list_1ary fn) [a]; 1720 1721/* Map a 2-ary function over a pair of objects. 1722 */ 1723map_binary fn a b = map_nary (list_2ary fn) [a, b]; 1724 1725/* Map a 3-ary function over three objects. 1726 */ 1727map_trinary fn a b c = map_nary (list_3ary fn) [a, b, c]; 1728 1729/* Map a 4-ary function over three objects. 1730 */ 1731map_quaternary fn a b c d = map_nary (list_4ary fn) [a, b, c, d]; 1732 1733/* Same as map_nary, but for lists. Handy for (eg.) implementing arith ops on 1734 * vectors and matricies. 1735 */ 1736map_nary_list fn args 1737 = fn args, lists == [] 1738 = map process [0, 1 .. shortest - 1] 1739{ 1740 // find all the list arguments 1741 lists = filter is_list args; 1742 1743 // what's the length of the shortest list arg? 1744 shortest = foldr1 min_pair (map len lists); 1745 1746 // process index n ... pull that member from each argument 1747 // recurse to handle application, so we work for nested lists too 1748 process n 1749 = map_nary_list fn (map (extract n) args) 1750 { 1751 extract n arg 1752 = arg?n, is_list arg 1753 = arg; 1754 } 1755} 1756 1757map_unaryl fn a = map_nary_list (list_1ary fn) [a]; 1758 1759map_binaryl fn a b = map_nary_list (list_2ary fn) [a, b]; 1760 1761/* Remove features smaller than x pixels across from an image. This used to be 1762 * rather complex ... convsep is now good enough to use. 1763 */ 1764smooth x image = convsep (matrix_gaussian_blur (to_real x * 2)) image; 1765 1766/* Chop up an image into a list of lists of smaller images. Pad edges with 1767 * black. 1768 */ 1769imagearray_chop tile_width tile_height hoverlap voverlap i 1770 = map chop' [0, vstep .. last_y] 1771{ 1772 width = get_width i; 1773 height = get_height i; 1774 bands = get_bands i; 1775 format = get_format i; 1776 type = get_type i; 1777 1778 tile_width' = to_real tile_width; 1779 tile_height' = to_real tile_height; 1780 hoverlap' = to_real hoverlap; 1781 voverlap' = to_real voverlap; 1782 1783 /* Unique pixels per tile. 1784 */ 1785 hstep = tile_width' - hoverlap'; 1786 vstep = tile_height' - voverlap'; 1787 1788 /* Number of tiles across and down. Remember the case where width == 1789 * hstep. 1790 */ 1791 across = (int) ((width - 1) / hstep) + 1; 1792 down = (int) ((height - 1) / vstep) + 1; 1793 1794 /* top/left of final tile. 1795 */ 1796 last_x = hstep * (across - 1); 1797 last_y = vstep * (down - 1); 1798 1799 /* How much do we need to pad by? 1800 */ 1801 sx = last_x + tile_width'; 1802 sy = last_y + tile_height'; 1803 1804 /* Expand image with black to pad size. 1805 */ 1806 pad = embed 0 0 0 sx sy i; 1807 1808 /* Chop up a row. 1809 */ 1810 chop' y 1811 = map chop'' [0, hstep .. last_x] 1812 { 1813 chop'' x = extract_area x y tile_width' tile_height' pad; 1814 } 1815} 1816 1817/* Reassemble image. 1818 */ 1819imagearray_assemble hoverlap voverlap il 1820 = (image_set_origin 0 0 @ foldl1 tbj @ map (foldl1 lrj)) il 1821{ 1822 lrj l r = insert (get_width l + hoverlap) 0 r l; 1823 tbj t b = insert 0 (get_height t + voverlap) b t; 1824} 1825 1826/* Generate an nxn identity matrix. 1827 */ 1828identity_matrix n 1829 = error "identity_matrix: n > 0", n < 1 1830 = map line [0 .. n - 1] 1831{ 1832 line p = take p [0, 0 ..] ++ [1] ++ take (n - p - 1) [0, 0 ..]; 1833} 1834 1835/* Incomplete gamma function Q(a, x) == 1 - P(a, x) 1836 1837 FIXME ... this is now a builtin, until we can get a nice List class 1838 (requires overloadable ':') 1839 1840gammq a x 1841 = error "bad args", x < 0 || a <= 0 1842 = 1 - gamser, x < a + 1 1843 = gammcf 1844{ 1845 gamser = (gser a x)?0; 1846 gammcf = (gcf a x)?0; 1847} 1848 */ 1849 1850/* Incomplete gamma function P(a, x) evaluated as series representation. Also 1851 * return ln(gamma(a)) ... nr in c, pp 218 1852 */ 1853gser a x 1854 = [gamser, gln] 1855{ 1856 gln = gammln a; 1857 gamser 1858 = error "bad args", x < 0 1859 = 0, x == 0 1860 = 1 // fix this 1861 { 1862 // maximum iterations 1863 maxit = 100; 1864 1865 ap = List [a + 1, a + 2 ...]; 1866 xoap = x / ap; 1867 1868 del = map product (prefixes xoap.value); 1869 1870 1871 1872 1873 1874 1875 1876/* 1877 del = map (multiply (1 / a)) (map product (prefixes xoap)) 1878 1879 del = 1880 1881 xap = iterate (multiply 1882 */ 1883 1884 /* Generate all prefixes of a list ... [1,2,3] -> [[1], [1, 2], [1, 2, 1885 * 3], [1, 2, 3, 4] ...] 1886 */ 1887 prefixes l = map (converse take l) [1..]; 1888 1889 } 1890} 1891 1892/* ln(gamma(xx)) ... nr in c, pp 214 1893 */ 1894gammln xx 1895 = gln 1896{ 1897 cof = [76.18009172947146, -86.50532032941677, 24.01409824083091, 1898 -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5]; 1899 y = take 6 (iterate (add 1) (xx + 1)); 1900 ser = 1.000000000190015 + sum (map2 divide cof y); 1901 tmp = xx + 0.5; 1902 tmp' = tmp - ((xx + 0.5) * log tmp); 1903 gln = -tmp + log (2.5066282746310005 * ser / xx); 1904} 1905 1906/* make a LUT from a scatter 1907 */ 1908buildlut x 1909 = Image (im_buildlut x), is_Matrix x && x.width > 1 1910 = im_buildlut (Matrix x), is_matrix x && is_list_len_more 1 x?0 1911 = error (_ "bad arguments to " ++ "buildlut"); 1912 1913/* Linear regression. Return a class with the stuff we need in. 1914 * from s15.2, p 665 NR in C 1915 */ 1916linreg xes yes 1917 = obj 1918{ 1919 obj = class { 1920 // in case we ever get shown in the workspace 1921 _vislevel = 2; 1922 1923 slope = sum [t * y :: [t, y] <- zip2 tes yes] / st2; 1924 intercept = (sy - sx * slope) / ss; 1925 1926 chi2 = sum [(y - intercept - slope * x) ** 2 :: [x, y] <- zip2 xes yes]; 1927 1928 siga = (chi2 / (ss - 2)) ** 0.5 * 1929 ((1 + sx ** 2 / (ss * st2)) / ss) ** 0.5; 1930 sigb = (chi2 / (ss - 2)) ** 0.5 * (1 / st2) ** 0.5; 1931 1932 // for compat with linregw, see below 1933 q = 1.0; 1934 } 1935 1936 ss = len xes; 1937 sx = sum xes; 1938 sy = sum yes; 1939 sxoss = sx / ss; 1940 1941 tes = [x - sxoss :: x <- xes]; 1942 st2 = sum [t ** 2 :: t <- tes]; 1943} 1944 1945/* Weighted linear regression. Xes, yes and a list of deviations. 1946 */ 1947linregw xes yes devs 1948 = obj 1949{ 1950 obj = class { 1951 // in case we ever get shown in the workspace 1952 _vislevel = 2; 1953 1954 slope = sum [(t * y) / sd :: [t, y, sd] <- zip3 tes yes devs] / st2; 1955 intercept = (sy - sx * slope) / ss; 1956 1957 chi2 = sum [((y - intercept - slope * x) / sd) ** 2 :: 1958 [x, y, sd] <- zip3 xes yes devs]; 1959 1960 siga = ((1 + sx * sx / (ss * st2)) / ss) ** 0.5; 1961 sigb = (1 / st2) ** 0.5; 1962 1963 q = gammq (0.5 * (len xes - 2)) (0.5 * chi2); 1964 } 1965 1966 wt = [sd ** -0.5 :: sd <- devs]; 1967 1968 ss = sum wt; 1969 sx = sum [x * w :: [x, w] <- zip2 xes wt]; 1970 sy = sum [y * w :: [y, w] <- zip2 yes wt]; 1971 sxoss = sx / ss; 1972 1973 tes = [(x - sxoss) / sd :: [x, sd] <- zip2 xes devs]; 1974 st2 = sum [t ** 2 :: t <- tes]; 1975} 1976 1977/* Clustering: pass in a list of points, repeatedly merge the 1978 * closest two points until no two points are closer than the threshold. 1979 * Return [merged-points, corresponding-weights]. A weight is a list of the 1980 * indexes we merged to make that point, ie. len weight == how significant 1981 * this point is. 1982 * 1983 * eg. 1984 * cluster 12 [152,154,155,42,159] == 1985 * [[155,42],[[1,2,0,4],[3]]] 1986 */ 1987cluster thresh points 1988 = oo_unary_function cluster_op points, is_class points 1989 // can't use [0..len points - 1], in case len points == 0 1990 = merge [points, map (converse cons []) (take (len points) [0 ..])], 1991 is_list points 1992 = error (_ "bad arguments to " ++ "cluster") 1993{ 1994 cluster_op = Operator "cluster" 1995 (cluster thresh) Operator_type.COMPOUND false; 1996 1997 merge x 1998 = x, m < 2 || d > thresh 1999 = merge [points', weights'] 2000 { 2001 [points, weights] = x; 2002 m = len points; 2003 2004 // generate indexes of all possible pairs, avoiding comparing a thing 2005 // to itself, and assuming that dist is reflexive 2006 // first index is always less than 2nd index 2007 // the +1,+2 makes sure we have an increasing generator, otherwise we 2008 // can get [3 .. 4] (for example), which will make a decreasing 2009 // sequence 2010 pairs = [[x, y] :: x <- [0 .. m - 1]; y <- [x + 1, x + 2 .. m - 1]]; 2011 2012 // distance function 2013 // arg is eg. [3,1], meaning get distance from point 3 to point 1 2014 dist x 2015 = abs (points?i - points?j) 2016 { 2017 [i, j] = x; 2018 } 2019 2020 // smallest distance, then the two points we merge 2021 p = minpos (map dist pairs); 2022 d = dist pairs?p; 2023 [i, j] = pairs?p; 2024 2025 // new point and new weight 2026 nw = weights?i ++ weights?j; 2027 np = (points?i * len weights?i + points?j * len weights?j) / len nw; 2028 2029 // remove element i from a list 2030 remove i l = take i l ++ drop (i + 1) l; 2031 2032 // remove two old points, add the new merged one 2033 // i < j (see "pairs", above) 2034 points' = np : remove i (remove j points); 2035 weights' = nw : remove i (remove j weights); 2036 } 2037} 2038 2039